* SUBROUTINE PA0GS1             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE:
* NUMERICAL COMPUTATION OF THE GRADIENT OF THE MODEL FUNCTION.
*
* PARAMETERS :
*  II  N  NUMBER OF VARIABLES.
*  II  KA  INDEF OF THE APPROXIMATED FUNCTION.
*  RI  X(N)  VECTOR OF VARIABLES.
*  RO  GA(N)  GRADIENT OF THE APPROXIMATED FUNCTION.
*  RI  FA  VALUE OF THE APPROXIMATED FUNCTION.
*  RI  ETA1  PRECISION OF THE COMPUTED VALUES.
*  IU  NAV  NUMBER OF APPROXIMATED FUNCTION EVALUATIONS.
*
      SUBROUTINE PA0GS1(N,KA,X,GA,FA,ETA1,NAV)
      INTEGER N,KA,NAV
      DOUBLE PRECISION X(*),GA(*),FA,ETA1
      DOUBLE PRECISION XSTEP,XTEMP,FTEMP,ETA
      INTEGER IVAR
      ETA=SQRT(ETA1)
      FTEMP=FA
      DO 4 IVAR=1,N
C
C     STEP SELECTION
C
      XSTEP=1.0D 0
      XSTEP=ETA*MAX(ABS(X(IVAR)),XSTEP)*SIGN(1.0D 0,X(IVAR))
      XTEMP=X(IVAR)
      X(IVAR)=X(IVAR)+XSTEP
      XSTEP=X(IVAR)-XTEMP
      NAV=NAV+1
      CALL FUN(N,KA,X,FA)
C
C     NUMERICAL DIFFERENTIATION
C
      GA(IVAR)=(FA-FTEMP)/XSTEP
      X(IVAR)=XTEMP
    4 CONTINUE
      FA=FTEMP
      RETURN
      END
* SUBROUTINE PA1SQ1             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION
* WHICH IS DEFINED AS A SUM OF SQUARES.
*
* PARAMETERS:
*  II  N  NUMBER OF VARIABLES.
*  RI  X(N)  VECTOR OF VARIABLES.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  AF(N)  VALUES OF THE APPROXIMATED FUNCTIONS.
*  RI  GA(NF)  GRADIENT OF THE APPROXIMATED FUNCTION.
*  RI  AG(N*N)  RECTANGULAR MATRIX WHICH IS USED FOR THE DIRECTION
*         VECTOR DETERMINATION.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  ETA1  PRECISION OF THE COMPUTES FUNCTION VALUES.
*  II  KDA  DEGREE OF COMPUTED DERIVATIVES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  IU  NFV  NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED.
*  IU  NFG  NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PA1SQ1(N,X,F,AF,GA,AG,G,ETA1,KDA,KD,LD,NFV,NFG)
      INTEGER N,KDA,KD,LD,NFV,NFG
      DOUBLE PRECISION X(*),F,AF(*),GA(*),AG(*),G(*),ETA1
      INTEGER KA,NAV
      DOUBLE PRECISION FA
      IF (KD.LE.LD) RETURN
      IF (KD.GE.0.AND.LD.LT.0) THEN
      F=0.0D0
      NFV=NFV+1
      ENDIF
      IF (KD.GE.1.AND.LD.LT.1) THEN
      CALL MXVSET(N,0.0D0,G)
      IF (KDA.GT.0) NFG=NFG+1
      ENDIF
      NAV=0
      DO 30 KA=1,N
      IF (KD.LT.0) GO TO 30
      IF (LD.GE.0) THEN
      FA = AF(KA)
      GO TO 10
      ELSE
      CALL FUN(N,KA,X,FA)
      AF(KA) = FA
      ENDIF
      IF (LD.GE.0) GO TO 10
      F=F+FA*FA
   10 IF (KD.LT.1) GO TO 30
      IF (KDA.GT.0) THEN
      CALL DFUN(N,KA,X,GA)
      ELSE
      CALL PA0GS1(N,KA,X,GA,FA,ETA1,NAV)
      ENDIF
      CALL MXVDIR(N,FA,GA,G,G)
      CALL MXVCOP(N,GA,AG((KA-1)*N+1))
   30 CONTINUE
      NFV=NFV+NAV/N
      IF (KD.GE.0.AND.LD.LT.0) F=0.5D0*F
      LD=KD
      RETURN
      END
* SUBROUTINE PA2SQ1             ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
*  COMPUTATION OF THE VALUE AND THE GRADIENT AND THE HESSIAN MATRIX
*  OF THE OBJECTIVE FUNCTION WHICH IS DEFINED AS A SUM OF SQUARES.
*
* PARAMETERS:
*  II  NF  NUMBER OF VARIABLES.
*  II  NA  NUMBER OF APPROXIMATED FUNCTIONS.
*  RI  GA(NF)  GRADIENT OF THE APPROXIMATED FUNCTION.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  HA(NF*(NF+1)/2)  HESSIAN MATRIX OF THE APPROXIMATED FUNCTION.
*  RO  H(NF*(NF+1)/2)  HESSIAN MATRIX OF THE OBJECTIVE FUNCTION.
*  RI  FA  VALUE OF THE SELECTED FUNCTION.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*
* COMMON DATA :
*  II  KDA  DEGREE OF ANALYTICALLY COMPUTED DERIVATIVES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  IU  NAV  NUMBER OF APPROXIMATED FUNCTION VALUES COMPUTED.
*  IU  NAG  NUMBER OF APPROXIMATED FUNCTION GRADIENTS COMPUTED.
*  IU  NAH  NUMBER OF APPROXIMATED FUNCTION HESSIAN MATRICES COMPUTED.
*  IU  NFV  NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED.
*  IU  NFG  NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED.
*  IU  NFH  NUMBER OF OBJECTIVE FUNCTION HESSIAN MATRICES COMPUTED.
*  IU  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*
* STATUS VARIABLES :
*  NS,ISP,TSS
*
* SUBPROGRAMS USED :
*  S  UYPRO1  PROLOGUE.
*  S  UYEPI1  EPILOGUE.
*  S  UYSET0  STATUS DEFINITION.
*  S  MXVSET  INITIATION OF A VECTOR.
*  S  MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S  MXDSMO  INITIATION OF A DENSE SYMMETRIC MATRIX.
*  S  MXDSMA  DENSE SYMMETRIC MATRIX AUGMENTED BY THE SCALED DENSE
*         SYMMETRIC MATRIX.
*  S  MXDSMU  CORRECTION OF A DENSE SYMMETRIC MATRIX.
*
      SUBROUTINE PA2SQ1(NF,NA,X,F,AF,GA,G,H,ETA1,KDA,KD,LD,NFV,NFG)
      INTEGER NF,NA,KDA,KD,LD,NFV,NFG
      DOUBLE PRECISION X(*),F,AF(*),GA(*),G(*),H(*),ETA1
      INTEGER KA,NAV
      DOUBLE PRECISION FA
      IF (KD.LE.LD) RETURN
      IF (KD.GE.0.AND.LD.LT.0) THEN
      F=0.0D0
      NFV=NFV+1
      ENDIF
      IF (KD.GE.1.AND.LD.LT.1) THEN
      CALL MXVSET(NF,0.0D0,G)
      IF (KDA.GT.0) NFG=NFG+1
      ENDIF
      IF (KD.GE.2.AND.LD.LT.2) CALL MXVSET(NF*(NF+1)/2,0.0D0,H)
      NAV=0
      DO 30 KA=1,NA
      IF (KD.LT.0) GO TO 30
      IF (LD.GE.0) THEN
      FA = AF(KA)
      GO TO 10
      ELSE
      CALL FUN(NF,KA,X,FA)
      AF(KA) = FA
      ENDIF
      IF (LD.GE.0) GO TO 10
      F=F+FA*FA
   10 IF (KD.LT.1) GO TO 30
      IF (KDA.GT.0) THEN
      CALL DFUN(NF,KA,X,GA)
      ELSE
      CALL PA0GS1(NF,KA,X,GA,FA,ETA1,NAV)
      ENDIF
      CALL MXVDIR(NF,FA,GA,G,G)
      IF (KD.LT.2) GO TO 30
      CALL MXDSMU(NF,H,1.0D0,GA)
   30 CONTINUE
      NFV=NFV+NAV/NA
      IF (KD.GE.0.AND.LD.LT.0) F=0.5D0*F
      LD=KD
      RETURN
      END
* SUBROUTINE PC1F01             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE CONSTRAINT FUNCTION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  FC  VALUE OF THE SELECTED CONSTRAINT FUNCTION.
*  RI  CF(NC)  VECTOR CONTAINING VALUES OF CONSTRAINT FUNCTIONS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  GC(NF)  GRADIENT OF THE SELECTED CONSTRAINT FUNCTION.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE GRADIENTS OF CONSTRAINT
*         FUNCTIONS.
*  RO  CMAX  MAXIMUM CONSTRAINT VIOLATION.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  II  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*
      SUBROUTINE PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD)
      DOUBLE PRECISION FC,CMAX
      INTEGER          KD,LD,NC,NF
      DOUBLE PRECISION CF(*),CG(*),CL(*),CU(*),GC(*),X(*)
      INTEGER          IC(*)
      DOUBLE PRECISION POM,TEMP
      INTEGER          KC
      IF (KD.LE.LD) RETURN
      IF (LD.LT.0) CMAX = 0.0D0
      DO 20 KC = 1,NC
          IF (KD.LT.0) GO TO 20
          IF (LD.GE.0) THEN
              FC = CF(KC)
              GO TO 10
          ELSE
              CALL CON(NF,KC,X,FC)
              CF(KC) = FC
          END IF
          IF (IC(KC).GT.0) THEN
              POM = 0.0D0
              TEMP = CF(KC)
              IF (IC(KC).EQ.1.OR.IC(KC).GE.3) POM = MIN(POM,
     +                                              TEMP - CL(KC))
              IF (IC(KC).EQ.2.OR.IC(KC).GE.3) POM = MIN(POM,
     +                                              CU(KC) - TEMP)
              IF (POM.LT.0.0D0) THEN
                  CMAX = MAX(CMAX,-POM)
              END IF
          END IF
   10     IF (KD.LT.1) GO TO 20
          IF (LD.GE.1) THEN
              CALL MXVCOP(NF,CG((KC - 1) * NF + 1),GC)
              GO TO 20
          ELSE
              CALL DCON(NF,KC,X,GC)
              CALL MXVCOP(NF,GC,CG((KC - 1) * NF + 1))
          END IF
   20 CONTINUE
      LD = KD
      RETURN
      END
* SUBROUTINE PF1F01                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
* PARAMETERS:
*  II  NF  NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  GF(NF)  GRADIENT OF THE MODEL FUNCTION.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  FF  VALUE OF THE MODEL FUNCTION.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  II  KD  DEGREE OF REQUIRED DERIVATIVES.
*  II  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  IEXT  TYPE OF EXTREMUM. IEXT=0-MINIMUM. IEXT=1-MAXIMUM.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*
      SUBROUTINE PF1F01(NF,X,GF,G,FF,F,KD,LD,IEXT)
      DOUBLE PRECISION F,FF
      INTEGER          IEXT,KD,LD,NF
      DOUBLE PRECISION GF(*),G(*),X(*)
      INTEGER          NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES
      COMMON           /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      IF (KD.LE.LD) RETURN
      IF (LD.GE.0) GO TO 10
      NFV = NFV + 1
      CALL OBJ(NF,X,FF)
      IF (IEXT.LE.0) THEN
          F =  FF
      ELSE
          F = -FF
      END IF
   10 IF (KD.LT.1) GO TO 20
      IF (LD.GE.1) GO TO 20
      NFG = NFG + 1
      CALL DOBJ(NF,X,GF)
      IF (IEXT.GT.0) THEN
          CALL MXVNEG(NF,GF,G)
      END IF
   20 LD = KD
      RETURN
      END
* SUBROUTINE PLADB0               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* NEW LINEAR CONSTRAINT OR A NEW SIMPLE BOUND IS ADDED TO THE
* ACTIVE SET.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  II  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   PLADR0  CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION
*         AFTER CONSTRAINT ADDITION.
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDRMV  COPY OF THE SELECTED COLUMN OF A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDRGR  PLANE ROTATION OF A TRANSPOSED DENSE RECTANGULAR MATRIX.
*  S   MXVORT  DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR
*         PLANE ROTATION.
*
      SUBROUTINE PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,
     & IER)
      INTEGER NF,N,ICA(*),INEW,NADD,IER
      DOUBLE PRECISION CG(*),CR(*),CZ(*),S(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION CK,CL
      INTEGER K,L,N1
      CALL PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      IF (IER.NE.0) RETURN
      IF (N.GT.0) THEN
      N1=N+1
      IF (INEW.GT.0) THEN
      CALL MXDRMM(NF,N1,CZ,CG((INEW-1)*NF+1),S)
      ELSE
      CALL MXDRMV(NF,N1,CZ,S,-INEW)
      ENDIF
      DO 1 L=1,N
      K=L+1
      CALL MXVORT(S(K),S(L),CK,CL,IER)
      CALL MXDRGR(NF,CZ,K,L,CK,CL,IER)
      IF (IER.LT.0) RETURN
    1 CONTINUE
      ENDIF
      IER=0
      RETURN
      END
* SUBROUTINE PLADB4               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* NEW LINEAR CONSTRAINT OR A NEW SIMPLE BOUND IS ADDED TO THE ACTIVE
* SET. TRANSFORMED HESSIAN MATRIX APPROXIMATION OR ITS INVERSION
* IS UPDATED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RU  H(NF*(NF+1)/2)  TRANSFORMED HESSIAN MATRIX APPROXIMATION OR
*         ITS INVERSION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  IU  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*         IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION.
*         IDECF=10-DIAGONAL MATRIX.
*  II  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   PLADR0  CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION
*         AFTER CONSTRAINT ADDITION.
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDRMV  COPY OF THE SELECTED COLUMN OF A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDRGR  PLANE ROTATION OF A TRANSPOSED DENSE RECTANGULAR MATRIX.
*         RECTANGULAR MATRIX.
*  S   MXDSMR  PLANE ROTATION OF A DENSE SYMMETRIC MATRIX.
*  S   MXVORT  DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR
*         PLANE ROTATION.
*
      SUBROUTINE PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,IDECF,
     & INEW,NADD,IER)
      INTEGER NF,N,ICA(*),IDECF,INEW,NADD,IER
      DOUBLE PRECISION CG(*),CR(*),CZ(*),H(*),S(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION CK,CL
      INTEGER I,J,K,L,N1
      IF (IDECF.NE.0.AND.IDECF.NE.9) THEN
      IER=-2
      RETURN
      ENDIF
      CALL PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      IF (IER.NE.0) RETURN
      IF (N.GT.0) THEN
      N1=N+1
      IF (INEW.GT.0) THEN
      CALL MXDRMM(NF,N1,CZ,CG((INEW-1)*NF+1),S)
      ELSE
      CALL MXDRMV(NF,N1,CZ,S,-INEW)
      ENDIF
      DO 1 L=1,N
      K=L+1
      CALL MXVORT(S(K),S(L),CK,CL,IER)
      CALL MXDRGR(NF,CZ,K,L,CK,CL,IER)
      CALL MXDSMR(N1,H,K,L,CK,CL,IER)
      IF (IER.LT.0) RETURN
    1 CONTINUE
      IF (IDECF.EQ.9) THEN
      L=N*(N+1)/2
      IF (H(L+N1).NE.0.0D 0) THEN
      CL=1.0D 0/H(L+N1)
      K=0
      DO 3 I=1,N
      CK=CL*H(L+I)
      DO 2 J=1,I
      K=K+1
      H(K)=H(K)-CK*H(L+J)
    2 CONTINUE
    3 CONTINUE
      ENDIF
      ENDIF
      ENDIF
      IER=0
      RETURN
      END
* SUBROUTINE PLADR0               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TRIANGULAR DECOMPOSITION OF KERNEL OF THE ORTHOGONAL PROJECTION
* IS UPDATED AFTER CONSTRAINT ADDITION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  II  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   MXSPRB  SPARSE BACK SUBSTITUTION.
*  S   MXVCOP  COPYING OF A VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      INTEGER NF,N,ICA(*),INEW,NADD,IER
      DOUBLE PRECISION CG(*),CR(*),S(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION MXVDOT
      INTEGER NCA,NCR,I,J,K,L
      IER=0
      IF (N.LE.0) IER=2
      IF (INEW.EQ.0) IER=3
      IF (IER.NE.0) RETURN
      NCA=NF-N
      NCR=NCA*(NCA+1)/2
      IF (INEW.GT.0) THEN
      CALL MXVCOP(NF,CG((INEW-1)*NF+1),S)
      GMAX=MXVDOT(NF,CG((INEW-1)*NF+1),S)
      DO 1 J=1,NCA
      L=ICA(J)
      IF (L.GT.0) THEN
      CR(NCR+J)=MXVDOT(NF,CG((L-1)*NF+1),S)
      ELSE
      I=-L
      CR(NCR+J)=S(I)
      ENDIF
    1 CONTINUE
      ELSE
      K=-INEW
      GMAX=1.0D 0
      DO 2 J=1,NCA
      L=ICA(J)
      IF (L.GT.0) THEN
      CR(NCR+J)=CG((L-1)*NF+K)*GMAX
      ELSE
      CR(NCR+J)=0.0D 0
      ENDIF
    2 CONTINUE
      ENDIF
      IF (NCA.EQ.0) THEN
      UMAX=GMAX
      ELSE
      CALL MXDPRB(NCA,CR,CR(NCR+1),1)
      UMAX=GMAX-MXVDOT(NCA,CR(NCR+1),CR(NCR+1))
      ENDIF
      IF (UMAX.LE.EPS7*GMAX) THEN
      IER=1
      RETURN
      ELSE
      N=N-1
      NCA=NCA+1
      NCR=NCR+NCA
      ICA(NCA)=INEW
      CR(NCR)=SQRT(UMAX)
      NADD=NADD+1
      ENDIF
      RETURN
      END
* SUBROUTINE PLLPB1             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE INITIAL FEASIBLE POINT AND THE LINEAR PROGRAMMING
* SUBROUTINE.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEAR CONSTRAINTS.
*  RU  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RU  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RA  CFD(NF)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  IO  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RO  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RO  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GO(NF)  SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  S(NF)  DIRECTION VECTOR.
*  II  MFP  TYPE OF FEASIBLE POINT. MFP=1-ARBITRARY FEASIBLE POINT.
*         MFP=2-OPTIMUM FEASIBLE POINT. MFP=3-REPEATED SOLUTION.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  RI  ETA9  MAXIMUM FOR REAL NUMBERS.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RI  EPS9  TOLERANCE FOR ACTIVITY OF CONSTRAINTS.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  IO  N  DIMENSION OF THE MANIFOLD DEFINED BY ACTIVE CONSTRAINTS.
*  IO  ITERL  TYPE OF FEASIBLE POINT. ITERL=1-ARBITRARY FEASIBLE POINT.
*         ITERL=2-OPTIMUM FEASIBLE POINT. ITERL=-1 FEASIBLE POINT DOES
*         NOT EXISTS. ITERL=-2 OPTIMUM FEASIBLE POINT DOES NOT EXISTS.
*
* SUBPROGRAMS USED :
*  S   PLINIT  DETERMINATION OF INITIAL POINT SATISFYING SIMPLE BOUNDS.
*  S   PLMAXL  MAXIMUM STEPSIZE USING LINEAR CONSTRAINTS.
*  S   PLMAXS  MAXIMUM STEPSIZE USING SIMPLE BOUNDS.
*  S   PLMAXT  MAXIMUM STEPSIZE USING TRUST REGION BOUNDS.
*  S   PLNEWL  IDENTIFICATION OF ACTIVE LINEAR CONSTRAINTS.
*  S   PLNEWS  IDENTIFICATION OF ACTIVE SIMPLE BOUNDS.
*  S   PLNEWT  IDENTIFICATION OF ACTIVE TRUST REGION BOUNDS.
*  S   PLDIRL  NEW VALUES OF CONSTRAINT FUNCTIONS.
*  S   PLDIRS  NEW VALUES OF VARIABLES.
*  S   PLSETC  INITIAL VALUES OF CONSTRAINT FUNCTIONS.
*  S   PLSETG  DETERMINATION OF THE FIRST PHASE GRADIENT VECTOR.
*  S   PLTRBG  DETERMINATION OF LAGRANGE MULTIPLIERS AND COMPUTATION
*  S   PLADB0  CONSTRAINT ADDITION.
*  S   PLRMB0  CONSTRAINT DELETION.
*  S   MXDCMM  PREMULTIPLICATION OF A VECTOR BY A DENSE RECTANGULAR
*         MATRIX STORED BY COLUMNS.
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A DENSE RECTANGULAR
*         MATRIX STORED BY ROWS.
*  S   MXDSMI  DETERMINATION OF THE INITIAL UNIT DENSE SYMMETRIC
*         MATRIX.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVINA  ABSOLUTE VALUES OF ELEMENTS OF AN INTEGER VECTOR.
*  S   MXVINC  UPDATE OF AN INTEGER VECTOR.
*  S   MXVIND  CHANGE OF THE INTEGER VECTOR FOR CONSTRAINT ADDITION.
*  S   MXVINT  CHANGE OF THE INTEGER VECTOR FOR TRUST REGION BOUND
*         ADDITION.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLLPB1(NF,NC,X,IX,XO,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,
     & CR,CZ,G,GO,S,MFP,KBF,KBC,ETA9,EPS7,EPS9,UMAX,GMAX,N,ITERL)
      INTEGER NF,NC,IX(*),IC(*),ICA(*),MFP,KBF,KBC,N,ITERL
      DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),CF(*),CFD(*),CL(*),
     & CU(*),CG(*),CR(*),CZ(*),G(*),GO(*),S(*),ETA9,EPS7,EPS9,
     & UMAX,GMAX
      DOUBLE PRECISION POM,CON,DMAX
      INTEGER NCA,NCR,NCZ,IPOM,I,K,IOLD,INEW,IER,KREM,KC,NRED
      INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      CON=ETA9
*
*     INITIATION
*
      CALL MXVCOP(NF,X,XO)
      IPOM=0
      NRED=0
      KREM=0
      ITERL=1
      DMAX=0.0D 0
      IF (MFP.EQ.3) GO TO 5
      IF (KBF.GT.0) CALL MXVINA(NF,IX)
*
*     SHIFT OF VARIABLES FOR SATISFYING SIMPLE BOUNDS
*
      CALL PLINIT(NF,X,IX,XL,XU,EPS9,KBF,INEW,ITERL)
      IF (ITERL.LT.0) THEN
      GO TO 6
      ENDIF
      N=0
      NCA=0
      NCZ=0
      DO 1 I=1,NF
      IF (KBF.GT.0.AND.IX(I).LT.0) THEN
      NCA=NCA+1
      ICA(NCA)=-I
      ELSE
      N=N+1
      CALL MXVSET(NF,0.0D 0,CZ(NCZ+1))
      CZ(NCZ+I)=1.0D 0
      NCZ=NCZ+NF
      ENDIF
    1 CONTINUE
      CALL MXDSMI(NCA,CR)
      IF (NC.GT.0) THEN
      CALL MXDRMM(NF,NC,CG,X,CF)
*
*     ADDITION OF ACTIVE CONSTRAINTS AND INITIAL CHECK OF FEASIBILITY
*
      CALL MXVINA(NC,IC)
C      IF (NF.GT.N) CALL PLSETC(NF,NC,X,XO,CF,IC,CG,S)
      DO 2 KC=1,NC
      IF (IC(KC).NE.0) THEN
      INEW=0
      CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      CALL MXVIND(IC,KC,IER)
      IF (IC(KC).LT.-10) IPOM=1
      ENDIF
    2 CONTINUE
      ENDIF
    3 IF (IPOM.EQ.1) THEN
*
*     CHECK OF FEASIBILITY AND UPDATE OF THE FIRST PHASE OBJECTIVE
*     FUNCTION
*
      CALL PLSETG(NF,NC,IC,CG,GO,INEW)
      IF (INEW.EQ.0) IPOM=0
      ENDIF
      IF (IPOM.EQ.0.AND.ITERL.EQ.0) THEN
*
*     FEASIBILITY ACHIEVED
*
      ITERL=1
      CALL MXVCOP(NF,G,GO)
      IF (MFP.EQ.1) GO TO 6
      ENDIF
*
*     LAGRANGE MULTIPLIERS AND REDUCED GRADIENT DETERMINATION
*
    5 CALL PLTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,GO,S,EPS7,GMAX,UMAX,IOLD)
      INEW=0
      IF (GMAX.EQ.0.0D 0) THEN
*
*     OPTIMUM ON A LINEAR MANIFOLD OBTAINED
*
      IF (IOLD.EQ.0) THEN
      IF (IPOM.EQ.0) THEN
*
*     OPTIMAL SOLUTION ACHIEVED
*
      ITERL= 2
      GO TO 6
      ELSE
      IPOM=0
      DO 22 KC=1,NC
      IF (IC(KC).LT.-10) THEN
      INEW=0
      CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      IF (IC(KC).LT.-10) IPOM=1
      ENDIF
   22 CONTINUE
      IF (IPOM.EQ.0) THEN
*
*     OPTIMAL SOLUTION ACHIEVED
*
      CALL MXVCOP(NF,GO,G)
      ITERL= 2
      GO TO 6
      ELSE
*
*     FEASIBLE SOLUTION DOES NOT EXIST
*
      CALL MXVCOP(NF,GO,G)
      ITERL=-1
      GO TO 6
      ENDIF
      ENDIF
      ELSE
*
*     CONSTRAINT DELETION
*
      CALL PLRMB0(NF,N,ICA,CG,CR,CZ,GO,S,IOLD,KREM,NREM,IER)
      KC=ICA(NF-N+1)
      IF (KC.GT.0) THEN
      IC(KC)=-IC(KC)
      ELSE
      K=-KC
      IX(K)=-IX(K)
      ENDIF
      DMAX=0.0D 0
      GO TO 5
      ENDIF
      ELSE
*
*     DIRECTION DETERMINATION
*
      NCA=NF-N
      NCR=NCA*(NCA+1)/2
      CALL MXDCMM(NF,N,CZ,S,CR(NCR+1))
      CALL MXVNEG(NF,CR(NCR+1),S)
*
*     STEPSIZE SELECTION
*
      POM=CON
      CALL PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,POM,KBC,KREM,INEW)
      CALL PLMAXS(NF,X,IX,XL,XU,S,POM,KBF,KREM,INEW)
      IF (INEW.EQ.0) THEN
      IF (IPOM.EQ.0) THEN
*
*     BOUNDED SOLUTION DOES NOT EXIST
*
      ITERL=-2
      ELSE
*
*     FEASIBLE SOLUTION DOES NOT EXIST
*
      ITERL=-3
      ENDIF
      GO TO 6
      ELSE
*
*     STEP REALIZATION
*
      CALL PLDIRS(NF,X,IX,S,POM,KBF)
      CALL PLDIRL(NC,CF,CFD,IC,POM,KBC)
*
*     CONSTRAINT ADDITION
*
      IF (INEW.GT.0) THEN
      KC=INEW
      INEW=0
      CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      CALL MXVIND(IC,KC,IER)
      ELSE IF (INEW+NF.GE.0) THEN
      I=-INEW
      INEW=0
      CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW)
      CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      CALL MXVIND(IX,I,IER)
      ENDIF
      DMAX=POM
      IF (DMAX.GT.0.0D 0) NRED=NRED+1
      GO TO 3
      ENDIF
      ENDIF
    6 CONTINUE
      RETURN
      END
* SUBROUTINE PLRMB0               ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* OLD LINEAR CONSTRAINT OR AN OLD SIMPLE BOUND IS REMOVED FROM THE
* ACTIVE SET.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  GN(NF)  TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION.
*  II  IOLD  INDEX OF THE OLD ACTIVE CONSTRAINT.
*  IO  KREM  AUXILIARY VARIABLE.
*  IU  NREM NUMBER OF CONSTRAINT DELETION.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   PLRMR0  CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION
*         AFTER CONSTRAINT DELETION.
*  S   MXDPRB  BACK SUBSTITUTION.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLRMB0(NF,N,ICA,CG,CR,CZ,G,GN,IOLD,KREM,NREM,IER)
      INTEGER NF,N,ICA(*),IOLD,KREM,NREM,IER
      DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*)
      DOUBLE PRECISION MXVDOT
      INTEGER NCA,NCR,NCZ,I,J,KC
      IER=0
      IF (N.EQ.NF) IER=2
      IF (IOLD.EQ.0) IER=3
      IF (IER.NE.0) RETURN
      NCA=NF-N
      NCR=NCA*(NCA-1)/2
      NCZ=N*NF
      CALL PLRMR0(NF,ICA,CR,CZ(NCZ+1),N,IOLD,KREM,IER)
      CALL MXVSET(NCA,0.0D 0,CZ(NCZ+1))
      CZ(NCZ+NCA)=1.0D 0
      CALL MXDPRB(NCA,CR,CZ(NCZ+1),-1)
      CALL MXVCOP(NCA,CZ(NCZ+1),CR(NCR+1))
      N=N+1
      CALL MXVSET(NF,0.0D 0,CZ(NCZ+1))
      DO 4 J=1,NCA
      KC=ICA(J)
      IF (KC.GT.0) THEN
      CALL MXVDIR(NF,CR(NCR+J),CG((KC-1)*NF+1),CZ(NCZ+1),CZ(NCZ+1))
      ELSE
      I=-KC
      CZ(NCZ+I)=CZ(NCZ+I)+CR(NCR+J)
      ENDIF
    4 CONTINUE
      GN(N)=MXVDOT(NF,CZ(NCZ+1),G)
      NREM=NREM+1
      IER=0
      RETURN
      END
* SUBROUTINE PLQDB1             ALL SYSTEMS                   97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  IO  N  DIMENSION OF THE MANIFOLD DEFINED BY ACTIVE CONSTRAINTS.
*  II  NC  NUMBER OF LINEAR CONSTRAINTS.
*  RU  X(NF)   VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IXA(NF)  VECTOR CONTAINING INFORMATION ON TRUST REGION ACTIVITY.
*  RI  XN(NF)  VECTOR OF SCALING FACTORS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RO  CFD(NC)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*            FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RO  CZ(NF)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RO  G(NF)  GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RO  GO(NF)  SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  H(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OR INVERSION OF THE
*         HESSIAN MATRIX APPROXIMATION.
*  RI  S(NF)  DIRECTION VECTOR.
*  RI  ETA2  TOLERANCE FOR POSITIVE DEFINITENESS OF THE HESSIAN MATRIX.
*  RI  ETA9  MAXIMUM FOR REAL NUMBERS.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RI  EPS9  TOLERANCE FOR ACTIVITY OF CONSTRAINTS.
*  RU  XDEL  TRUST REGION BOUND.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  II  MFP  TYPE OF FEASIBLE POINT. MFP=1-ARBITRARY FEASIBLE POINT.
*         MFP=2-OPTIMUM FEASIBLE POINT. MFP=3-REPEATED SOLUTION.
*
* COMMON DATA :
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  NORMF  SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED.
*         NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY.
*         NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER.
*  IU  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*         IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION.
*         IDECF=10-DIAGONAL MATRIX.
*  IU  NDECF  NUMBER OF DECOMPOSITIONS.
*  IO  ITERQ  TYPE OF FEASIBLE POINT. ITERQ=1-ARBITRARY FEASIBLE POINT.
*         ITERQ=2-OPTIMUM FEASIBLE POINT. ITERQ=-1 FEASIBLE POINT DOES
*         NOT EXISTS. ITERQ=-2 OPTIMUM FEASIBLE POINT DOES NOT EXISTS.
*
* SUBPROGRAMS USED :
*  S   PLMINS  DETERMINATION OF THE NEW ACTIVE SIMPLE BOUND.
*  S   PLMINL  DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT.
*  S   PLMINT  DETERMINATION OF THE NEW ACTIVE TRUST REGION BOUND.
*  S   PLADR1  ADDITION OF A NEW ACTIVE CONSTRAINT.
*  S   PLRMR0  CONSTRAIN DELETION.
*  S   PLSOB1  TRANSFORMATION OF THE LOCAL SOLUTION TO THE SOLUTION
*         OF THE ORIGINAL QP PROBLEM.
*  S   MXDPGF  GILL-MURRAY DECOMPOSITION OF A DENSE SYMMETRIC MATRIX.
*  S   MXDPGB  BACK SUBSTITUTION AFTER GILL-MURRAY DECOMPOSITION.
*  S   MXDPRB  BACK SUBSTITUTION.
*  S   MXDSMM  MATRIX VECTOR PRODUCT.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVINA  ABSOLUTE VALUES OF ELEMENTS OF AN INTEGER VECTOR.
*  S   MXVINV  CHANGE OF AN INTEGER VECTOR AFTER CONSTRAINT ADDITION.
*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*
      SUBROUTINE PLQDB1(NF,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,CR,CZ,
     & G,GO,H,S,MFP,KBF,KBC,IDECF,ETA2,ETA9,EPS7,EPS9,UMAX,GMAX,N,ITERQ)
      INTEGER NF,NC,IX(*),IC(*),ICA(*),MFP,KBF,KBC,IDECF,N,ITERQ
      DOUBLE PRECISION X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*),CG(*),
     & CR(*),CZ(*),G(*),GO(*),H(*),S(*),ETA2,ETA9,EPS7,EPS9,UMAX,GMAX
      DOUBLE PRECISION CON,TEMP,STEP,STEP1,STEP2,DMAX,PAR,SNORM
      INTEGER NCA,NCR,I,J,K,IOLD,JOLD,INEW,JNEW,KNEW,INF,IER,KREM,KC,
     & NRED
      INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      CON=ETA9
      IF (IDECF.LT.0) IDECF=1
      IF (IDECF.EQ.0) THEN
C
C     GILL-MURRAY DECOMPOSITION
C
      TEMP=ETA2
      CALL MXDPGF(NF,H,INF,TEMP,STEP)
      NDEC=NDEC+1
      IDECF=1
      ENDIF
      IF (IDECF.GE.2.AND.IDECF.LE.8) THEN
      ITERQ=-10
      RETURN
      ENDIF
C
C     INITIATION
C
      NRED=0
      JOLD=0
      JNEW=0
      ITERQ=0
      DMAX=0.0D0
      IF (MFP.EQ.3) GO TO 1
      N=NF
      NCA=0
      NCR=0
      IF (KBF.GT.0) CALL MXVINA(NF,IX)
      IF (KBC.GT.0) CALL MXVINA(NC,IC)
C
C     DIRECTION DETERMINATION
C
    1 CALL MXVNEG(NF,GO,S)
      DO 2 J=1,NCA
      KC=ICA(J)
      IF (KC.GT.0) THEN
      CALL MXVDIR(NF,CZ(J),CG((KC-1)*NF+1),S,S)
      ELSE
      K=-KC
      S(K)=S(K)+CZ(J)
      ENDIF
    2 CONTINUE
      CALL MXVCOP(NF,S,G)
      IF (IDECF.EQ.1) THEN
      CALL MXDPGB(NF,H,S,0)
      ELSE
      CALL MXDSMM(NF,H,G,S)
      ENDIF
      IF (ITERQ.EQ.3) GO TO 7
C
C     CHECK OF FEASIBILITY
C
      INEW=0
      PAR=0.0D0
      CALL PLMINN(NF,NC,CF,CFD,IC,CL,CU,CG,S,EPS9,PAR,KBC,INEW,KNEW)
      CALL PLMINS(NF,IX,X,XL,XU,S,KBF,INEW,KNEW,EPS9,PAR)
      IF (INEW.EQ.0) THEN
C
C     SOLUTION ACHIEVED
C
      CALL MXVNEG(NF,G,G)
      ITERQ=2
      GO TO 7
      ELSE
      SNORM=0.0D0
      ENDIF
    4 IER=0
C
C     STEPSIZE DETERMINATION
C
      CALL PLADR1(NF,N,ICA,CG,CR,H,S,G,EPS7,GMAX,UMAX,IDECF,INEW,
     & NADD,IER,1)
      CALL MXDPRB(NCA,CR,G,-1)
      IF (KNEW.LT.0) CALL MXVNEG(NCA,G,G)
C
C     PRIMAL STEPSIZE
C
      IF (IER.NE.0) THEN
      STEP1=CON
      ELSE
      STEP1=-PAR/UMAX
      ENDIF
C
C     DUAL STEPSIZE
C
      IOLD=0
      STEP2=CON
      DO 5 J=1,NCA
      KC=ICA(J)
      IF (KC.GE.0) THEN
      K=IC(KC)
      ELSE
      I=-KC
      K=IX(I)
      ENDIF
      IF (K.LE.-5) THEN
      ELSE IF ((K.EQ.-1.OR.K.EQ.-3.).AND.G(J).LE.0.0D0) THEN
      ELSE IF ((K.EQ.-2.OR.K.EQ.-4.).AND.G(J).GE.0.0D0) THEN
      ELSE
      TEMP=CZ(J)/G(J)
      IF (STEP2.GT.TEMP) THEN
      IOLD=J
      STEP2=TEMP
      ENDIF
      ENDIF
    5 CONTINUE
C
C     FINAL STEPSIZE
C
      STEP=MIN(STEP1,STEP2)
      IF (STEP.GE.CON) THEN
C
C     FEASIBLE SOLUTION DOES NOT EXIST
C
      ITERQ=-1
      GO TO 7
      ENDIF
C
C     NEW LAGRANGE MULTIPLIERS
C
      DMAX=STEP
      CALL MXVDIR(NCA,-STEP,G,CZ,CZ)
      SNORM=SNORM+SIGN(1,KNEW)*STEP
      PAR=PAR-(STEP/STEP1)*PAR
      IF (STEP.EQ.STEP1) THEN
      IF (N.LE.0) THEN
C
C     IMPOSSIBLE SITUATION
C
      ITERQ=-5
      GO TO 7
      ENDIF
C
C     CONSTRAINT ADDITION
C
      IF (IER.EQ.0) THEN
      N=N-1
      NCA=NCA+1
      NCR=NCR+NCA
      CZ(NCA)=SNORM
      ENDIF
      IF (INEW.GT.0) THEN
      KC=INEW
      CALL MXVINV(IC,KC,KNEW)
      ELSE IF (ABS(KNEW).EQ.1) THEN
      I=-INEW
      CALL MXVINV(IX,I,KNEW)
      ELSE
      I=-INEW
      IF (KNEW.GT.0) IX(I)=-3
      IF (KNEW.LT.0) IX(I)=-4
      ENDIF
      NRED=NRED+1
      NADD=NADD+1
      JNEW=INEW
      JOLD=0
      GO TO 1
      ELSE
C
C     CONSTRAINT DELETION
C
      DO 6 J=IOLD,NCA-1
      CZ(J)=CZ(J+1)
    6 CONTINUE
      CALL PLRMF0(NF,NC,IX,IC,ICA,CR,IC,G,N,IOLD,KREM,IER)
      NCR=NCR-NCA
      NCA=NCA-1
      JOLD=IOLD
      JNEW=0
      IF (KBC.GT.0) CALL MXVINA(NC,IC)
      IF (KBF.GT.0) CALL MXVINA(NF,IX)
      DO 8 J=1,NCA
      KC=ICA(J)
      IF (KC.GT.0) THEN
      IC(KC)=-IC(KC)
      ELSE
      KC=-KC
      IX(KC)=-IX(KC)
      ENDIF
    8 CONTINUE
      GO TO 4
      ENDIF
    7 CONTINUE
      RETURN
      END
* SUBROUTINE PLADR1               ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TRIANGULAR DECOMPOSITION OF KERNEL OF THE GENERAL PROJECTION
* IS UPDATED AFTER CONSTRAINT ADDITION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  H(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OR INVERSION OF THE
*         HESSIAN MATRIX APPROXIMATION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RO  G(NF)  VECTOR USED IN THE DUAL RANGE SPACE QUADRATIC PROGRAMMING
*         METHOD.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  RO  E  AUXILIARY VARIABLE.
*  RI  T  AUXILIARY VARIABLE.
*  IU  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*         IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION.
*         IDECF=10-DIAGONAL MATRIX.
*  II  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  IER  ERROR INDICATOR.
*  II  JOB  SPECIFICATION OF COMPUTATION. OUTPUT VECTOR G IS NOT OR IS
*         COMPUTED IN CASE WHEN N.LE.0 IF JOB=0 OR JOB=1 RESPECTIVELY.
*
* SUBPROGRAMS USED :
*  S   MXDPGB  BACK SUBSTITUTION.
*  S   MXDPRB  BACK SUBSTITUTION.
*  S   MXDSMM  MATRIX-VECTOR PRODUCT.
*  S   MXDSMV  COPYING OF A ROW OF DENSE SYMMETRIC MATRIX.
*  S   MXVCOP  COPYING OF A VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLADR1(NF,N,ICA,CG,CR,H,S,G,EPS7,GMAX,UMAX,IDECF,
     & INEW,NADD,IER,JOB)
      INTEGER NF,N,ICA(*),IDECF,INEW,NADD,IER,JOB
      DOUBLE PRECISION CG(*),CR(*),H(*),S(*),G(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION MXVDOT
      INTEGER NCA,NCR,JCG,J,K,L
      IER=0
      IF (JOB.EQ.0.AND.N.LE.0) IER=2
      IF (INEW.EQ.0) IER=3
      IF (IDECF.NE.1.AND.IDECF.NE.9) IER=-2
      IF (IER.NE.0) RETURN
      NCA=NF-N
      NCR=NCA*(NCA+1)/2
      IF (INEW.GT.0) THEN
      JCG=(INEW-1)*NF+1
      IF (IDECF.EQ.1) THEN
      CALL MXVCOP(NF,CG(JCG),S)
      CALL MXDPGB(NF,H,S,0)
      ELSE
      CALL MXDSMM(NF,H,CG(JCG),S)
      ENDIF
      GMAX=MXVDOT(NF,CG(JCG),S)
      ELSE
      K=-INEW
      IF (IDECF.EQ.1) THEN
      CALL MXVSET(NF,0.0D0,S)
      S(K)=1.0D 0
      CALL MXDPGB(NF,H,S,0)
      ELSE
      CALL MXDSMV(NF,H,S,K)
      ENDIF
      GMAX=S(K)
      ENDIF
      DO 1 J=1,NCA
      L=ICA(J)
      IF (L.GT.0) THEN
      G(J)=MXVDOT(NF,CG((L-1)*NF+1),S)
      ELSE
      L=-L
      G(J)=S(L)
      ENDIF
    1 CONTINUE
      IF (N.EQ.0) THEN
      CALL MXDPRB(NCA,CR,G,1)
      UMAX=0.0D0
      IER=2
      RETURN
      ELSE IF (NCA.EQ.0) THEN
      UMAX=GMAX
      ELSE
      CALL MXDPRB(NCA,CR,G,1)
      UMAX=GMAX-MXVDOT(NCA,G,G)
      CALL MXVCOP(NCA,G,CR(NCR+1))
      ENDIF
      IF (UMAX.LE.EPS7*GMAX) THEN
      IER=1
      RETURN
      ELSE
      NCA=NCA+1
      NCR=NCR+NCA
      ICA(NCA)=INEW
      CR(NCR)=SQRT(UMAX)
      IF (JOB.EQ.0) THEN
      N=N-1
      NADD=NADD+1
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PLDIRL               ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE NEW VALUES OF THE CONSTRAINT FUNCTIONS.
*
* PARAMETERS :
*  II  NC  NUMBER OF CONSTRAINTS.
*  RU  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RI  CFD(NF)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  STEP  CURRENT STEPSIZE.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*
      SUBROUTINE PLDIRL(NC,CF,CFD,IC,STEP,KBC)
      INTEGER NC,IC(*),KBC
      DOUBLE PRECISION CF(*),CFD(*),STEP
      INTEGER KC
      IF (KBC.GT.0) THEN
      DO 1 KC=1,NC
      IF (IC(KC).GE.0.AND.IC(KC).LE.10) THEN
      CF(KC)=CF(KC)+STEP*CFD(KC)
      ELSE IF (IC(KC).LT.-10) THEN
      CF(KC)=CF(KC)+STEP*CFD(KC)
      ENDIF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLDIRS               ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE NEW VECTOR OF VARIABLES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RU  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  S(NF)  DIRECTION VECTOR.
*  RI  STEP  CURRENT STEPSIZE.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*
      SUBROUTINE PLDIRS(NF,X,IX,S,STEP,KBF)
      INTEGER NF,IX(*),KBF
      DOUBLE PRECISION X(*),S(*),STEP
      INTEGER I
      DO 1 I=1,NF
      IF (KBF.LE.0) THEN
      X(I)=X(I)+STEP*S(I)
      ELSE IF (IX(I).GE.0.AND.IX(I).LE.10) THEN
      X(I)=X(I)+STEP*S(I)
      ELSE IF (IX(I).LT.-10) THEN
      X(I)=X(I)+STEP*S(I)
      ENDIF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PLINIT             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE INITIAL POINT WHICH SATISFIES SIMPLE BOUNDS.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RU  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IO  IND  INDICATOR. IF IND.NE.0 THEN TRUST REGION BOUNDS CANNOT
*         BE SATISFIED.
*
* SUBPROGRAMS USED :
*  S   PLNEWS  TEST ON ACTIVITY OF A GIVEN SIMPLE BOUND.
*
      SUBROUTINE PLINIT(NF,X,IX,XL,XU,EPS9,KBF,INEW,IND)
      INTEGER NF,IX(*),KBF,INEW,IND
      DOUBLE PRECISION X(*),XL(*),XU(*),EPS9
      INTEGER I
      IND=0
      IF (KBF.GT.0) THEN
      DO 1 I=1,NF
      CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW)
      IF (IX(I).LT.5) THEN
      ELSE IF (IX(I).EQ.5) THEN
      IX(I)=-5
      ELSE IF (IX(I).EQ.11.OR.IX(I).EQ.13) THEN
      X(I)=XL(I)
      IX(I)=10-IX(I)
      ELSE IF (IX(I).EQ.12.OR.IX(I).EQ.14) THEN
      X(I)=XU(I)
      IX(I)=10-IX(I)
      ENDIF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLMAXL               ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE MAXIMUM STEPSIZE USING LINEAR CONSTRAINTS.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CURRENT LINEAR CONSTRAINTS.
*  RI  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS.
*  RO  CFD(NF)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  S(NF)  DIRECTION VECTOR.
*  RO  STEP  MAXIMUM STEPSIZE.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  KREM  INDICATION OF LINEARLY DEPENDENT GRADIENTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE FUNCTION.
*
* SUBPROGRAMS USED :
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,STEP,KBC,KREM,INEW)
      INTEGER NF,NC,IC(*),KBC,KREM,INEW
      DOUBLE PRECISION CF(*),CFD(*),CL(*),CU(*),CG(*),S(*),STEP
      DOUBLE PRECISION TEMP,MXVDOT
      INTEGER JCG,KC
      IF (KBC.GT.0) THEN
      JCG=1
      DO 1 KC=1,NC
      IF (KREM.GT.0.AND.IC(KC).GT.10) IC(KC)=IC(KC)-10
      IF (IC(KC).GT.0.AND.IC(KC).LE.10) THEN
      TEMP=MXVDOT(NF,CG(JCG),S)
      CFD(KC)=TEMP
      IF (TEMP.LT.0.0D 0) THEN
      IF (IC(KC).EQ.1.OR.IC(KC).GE.3) THEN
      TEMP=(CL(KC)-CF(KC))/TEMP
      IF (TEMP.LE.STEP) THEN
      INEW=KC
      STEP=TEMP
      ENDIF
      ENDIF
      ELSE IF (TEMP.GT.0.0D 0) THEN
      IF (IC(KC).EQ.2.OR.IC(KC).GE.3) THEN
      TEMP=(CU(KC)-CF(KC))/TEMP
      IF (TEMP.LE.STEP) THEN
      INEW=KC
      STEP=TEMP
      ENDIF
      ENDIF
      ENDIF
      ELSE IF (IC(KC).LT.-10) THEN
      TEMP=MXVDOT(NF,CG(JCG),S)
      CFD(KC)=TEMP
      IF (TEMP.GT.0.0D 0) THEN
      IF (IC(KC).EQ.-11.OR.IC(KC).EQ.-13.OR.IC(KC).EQ.-15) THEN
      TEMP=(CL(KC)-CF(KC))/TEMP
      IF (TEMP.LE.STEP) THEN
      INEW=KC
      STEP=TEMP
      ENDIF
      ENDIF
      ELSE IF (TEMP.LT.0.0D 0) THEN
      IF (IC(KC).EQ.-12.OR.IC(KC).EQ.-14.OR.IC(KC).EQ.-16) THEN
      TEMP=(CU(KC)-CF(KC))/TEMP
      IF (TEMP.LE.STEP) THEN
      INEW=KC
      STEP=TEMP
      ENDIF
      ENDIF
      ENDIF
      ENDIF
      JCG=JCG+NF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLMAXS               ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE MAXIMUM STEPSIZE USING THE SIMPLE BOUNDS
* FOR VARIABLES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  S(NF)  DIRECTION VECTOR.
*  RO  STEP  MAXIMUM STEPSIZE.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  IO  KREM  INDICATION OF LINEARLY DEPENDENT GRADIENTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*
      SUBROUTINE PLMAXS(NF,X,IX,XL,XU,S,STEP,KBF,KREM,INEW)
      INTEGER NF,IX(*),KBF,KREM,INEW
      DOUBLE PRECISION X(*),XL(*),XU(*),S(*),STEP
      DOUBLE PRECISION TEMP
      INTEGER I
      IF (KBF.GT.0) THEN
      DO 1 I=1,NF
      IF (KREM.GT.0.AND.IX(I).GT.10) IX(I)=IX(I)-10
      IF (IX(I).GT.0.AND.IX(I).LE.10) THEN
      IF (S(I).LT.0.0D 0) THEN
      IF (IX(I).EQ.1.OR.IX(I).GE.3) THEN
      TEMP=(XL(I)-X(I))/S(I)
      IF (TEMP.LE.STEP) THEN
      INEW=-I
      STEP=TEMP
      ENDIF
      ENDIF
      ELSE IF (S(I).GT.0.0D 0) THEN
      IF (IX(I).EQ.2.OR.IX(I).GE.3) THEN
      TEMP=(XU(I)-X(I))/S(I)
      IF (TEMP.LE.STEP) THEN
      INEW=-I
      STEP=TEMP
      ENDIF
      ENDIF
      ENDIF
      ENDIF
    1 CONTINUE
      ENDIF
      KREM=0
      RETURN
      END
* SUBROUTINE PLNEWL             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TEST ON ACTIVITY OF A GIVEN LINEAR CONSTRAINT.
*
* PARAMETERS :
*  II  KC  INDEX OF A GIVEN CONSTRAINT.
*  RI  CF(NC)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  IU  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*
      SUBROUTINE PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      INTEGER KC,IC(*),INEW
      DOUBLE PRECISION CF(*),CL(*),CU(*),EPS9
      DOUBLE PRECISION TEMP
      IF (IC(KC).LT.-10) IC(KC)=-IC(KC)-10
      IF (IC(KC).LE.0) THEN
      ELSE IF (IC(KC).EQ.1) THEN
      TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0)
      IF (CF(KC).GT.CL(KC)+TEMP) THEN
      ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN
      IC(KC)=11
      INEW=KC
      ELSE
      IC(KC)=-11
      ENDIF
      ELSE IF (IC(KC).EQ.2) THEN
      TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0)
      IF (CF(KC).LT.CU(KC)-TEMP) THEN
      ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN
      IC(KC)=12
      INEW=KC
      ELSE
      IC(KC)=-12
      ENDIF
      ELSE IF (IC(KC).EQ.3.OR.IC(KC).EQ.4) THEN
      TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0)
      IF (CF(KC).GT.CL(KC)+TEMP) THEN
      TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0)
      IF (CF(KC).LT.CU(KC)-TEMP) THEN
      ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN
      IC(KC)=14
      INEW=KC
      ELSE
      IC(KC)=-14
      ENDIF
      ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN
      IC(KC)=13
      INEW=KC
      ELSE
      IC(KC)=-13
      ENDIF
      ELSE IF (IC(KC).EQ.5.OR.IC(KC).EQ.6) THEN
      TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0)
      IF (CF(KC).GT.CL(KC)+TEMP) THEN
      TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0)
      IF (CF(KC).LT.CU(KC)-TEMP) THEN
      ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN
      IC(KC)=16
      INEW=KC
      ELSE
      IC(KC)=-16
      ENDIF
      ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN
      IC(KC)=15
      INEW=KC
      ELSE
      IC(KC)=-15
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PLMINN             ALL SYSTEMS                   97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  CF(NC)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RO  CFD(NC)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*            FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  S(NF)  DIRECTION VECTOR.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  RA  PAR  AUXILIARY VARIABLE.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IO  KNEW  SIGNUM OF THE NEW ACTIVE NORMAL.
*
* SUBPROGRAMS USED :
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLMINN(NF,NC,CF,CFD,IC,CL,CU,CG,S,EPS9,PAR,KBC,INEW,
     & KNEW)
      INTEGER NF,NC,IC(*),KBC,INEW,KNEW
      DOUBLE PRECISION CF(*),CFD(*),CL(*),CU(*),CG(*),S(*),EPS9,PAR
      DOUBLE PRECISION TEMP,POM,MXVDOT
      INTEGER JCG,KC
      IF (KBC.GT.0) THEN
      JCG=1
      DO 1 KC=1,NC
      IF (IC(KC).GT.0) THEN
      TEMP=MXVDOT(NF,CG(JCG),S)
      CFD(KC)=TEMP
      TEMP=CF(KC)+TEMP
      IF (IC(KC).EQ.1.OR.IC(KC).GE.3) THEN
      POM=TEMP-CL(KC)
      IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(CL(KC)),1.0D 0))) THEN
      INEW=KC
      KNEW= 1
      PAR=POM
      ENDIF
      ENDIF
      IF (IC(KC).EQ.2.OR.IC(KC).GE.3) THEN
      POM=CU(KC)-TEMP
      IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(CU(KC)),1.0D 0))) THEN
      INEW=KC
      KNEW=-1
      PAR=POM
      ENDIF
      ENDIF
      ENDIF
      JCG=JCG+NF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLMINS             ALL SYSTEMS                   91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE NEW ACTIVE SIMPLE BOUND.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XO(NF)  SAVED VECTOR OF VARIABLES.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  S(NF)  DIRECTION VECTOR.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IO  KNEW  SIGNUM OF THE NEW NORMAL.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  RA  PAR  AUXILIARY VARIABLE.
*
      SUBROUTINE PLMINS(NF,IX,XO,XL,XU,S,KBF,INEW,KNEW,EPS9,PAR)
      DOUBLE PRECISION EPS9,PAR
      INTEGER          INEW,KBF,KNEW,NF
      DOUBLE PRECISION S(*),XL(*),XO(*),XU(*)
      INTEGER          IX(*)
      DOUBLE PRECISION POM,TEMP
      INTEGER          I
      IF (KBF.GT.0) THEN
          DO 10 I = 1,NF
              IF (IX(I).GT.0) THEN
                  TEMP = 1.0D0
                  IF (IX(I).EQ.1 .OR. IX(I).GE.3) THEN
                      POM = XO(I) + S(I)*TEMP - XL(I)
                      IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(XL(I)),
     +                    TEMP))) THEN
                          INEW = -I
                          KNEW = 1
                          PAR = POM
                      END IF
                  END IF
                  IF (IX(I).EQ.2 .OR. IX(I).GE.3) THEN
                      POM = XU(I) - S(I)*TEMP - XO(I)
                      IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(XU(I)),
     +                    TEMP))) THEN
                          INEW = -I
                          KNEW = -1
                          PAR = POM
                      END IF
                  END IF
              END IF
   10     CONTINUE
      END IF
      RETURN
      END
* SUBROUTINE PLNEWS             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TEST ON ACTIVITY OF A GIVEN SIMPLE BOUND.
*
* PARAMETERS :
*  RI  X(NF)  VECTOR OF VARIABLES.
*  IU  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  II  I  INDEX OF TESTED SIMPLE BOUND.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*
      SUBROUTINE PLNEWS(X,IX,XL,XU,EPS9,I,INEW)
      INTEGER IX(*),I,INEW
      DOUBLE PRECISION X(*),XL(*),XU(*),EPS9
      DOUBLE PRECISION TEMP
      TEMP=1.0D 0
      IF (IX(I).LE.0) THEN
      ELSE IF (IX(I).EQ.1) THEN
      IF (X(I).LE.XL(I)+EPS9*MAX(ABS(XL(I)),TEMP)) THEN
      IX(I)=11
      INEW=-I
      ENDIF
      ELSE IF (IX(I).EQ.2) THEN
      IF (X(I).GE.XU(I)-EPS9*MAX(ABS(XU(I)),TEMP)) THEN
      IX(I)=12
      INEW=-I
      ENDIF
      ELSE IF (IX(I).EQ.3.OR.IX(I).EQ.4) THEN
      IF (X(I).LE.XL(I)+EPS9*MAX(ABS(XL(I)),TEMP)) THEN
      IX(I)=13
      INEW=-I
      ENDIF
      IF (X(I).GE.XU(I)-EPS9*MAX(ABS(XU(I)),TEMP)) THEN
      IX(I)=14
      INEW=-I
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PLREDL               ALL SYSTEMS                   98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TRANSFORMATION OF THE INCOMPATIBLE QUADRATIC PROGRAMMING SUBPROBLEM.
*
* PARAMETERS :
*  II  NC  NUMBER OF CURRENT LINEAR CONSTRAINTS.
*  RI  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*
      SUBROUTINE PLREDL(NC,CF,IC,CL,CU,KBC)
      INTEGER NC,IC(NC),KBC,K
      DOUBLE PRECISION CF(*),CL(*),CU(*)
      DOUBLE PRECISION TEMP
      INTEGER KC
      IF (KBC.GT.0) THEN
      DO 1 KC=1,NC
      K=IC(KC)
      IF (ABS(K).EQ.1.OR.ABS(K).EQ.3.OR.ABS(K).EQ.4) THEN
      TEMP=(CF(KC)-CL(KC))
      IF (TEMP.LT.0) CF(KC)=CL(KC)+0.1D 0*TEMP
      ENDIF
      IF (ABS(K).EQ.2.OR.ABS(K).EQ.3.OR.ABS(K).EQ.4) THEN
      TEMP=(CF(KC)-CU(KC))
      IF (TEMP.GT.0) CF(KC)=CU(KC)+0.1D 0*TEMP
      ENDIF
      IF (ABS(K).EQ.5.OR.ABS(K).EQ.6) THEN
      TEMP=(CF(KC)-CL(KC))
      CF(KC)=CL(KC)+0.1D 0*TEMP
      ENDIF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLRMF0             ALL SYSTEMS                   91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* OPERATIONS AFTER CONSTRAINT DELETION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IA(NA)  VECTOR CONTAINING TYPES OF DEVIATIONS.
*  IU  IAA(NF+1)  VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS.
*  RU  AR((NF+1)*(NF+2)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RA  S(NF+1)  AUXILIARY VECTOR.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  IOLD  INDEX OF THE OLD ACTIVE CONSTRAINT.
*  IO  KREM  AUXILIARY VARIABLE.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   PLRMR0  CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION
*         AFTER CONSTRAINT DELETION.
*
      SUBROUTINE PLRMF0(NF,NC,IX,IA,IAA,AR,IC,S,N,IOLD,KREM,IER)
      INTEGER          IER,IOLD,KREM,N,NC,NF
      DOUBLE PRECISION AR(*),S(*)
      INTEGER          IA(*),IAA(*),IC(*),IX(*)
      INTEGER          NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES
      INTEGER          L
      COMMON           /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      CALL PLRMR0(NF,IAA,AR,S,N,IOLD,KREM,IER)
      N = N + 1
      NREM = NREM + 1
      L = IAA(NF-N+1)
      IF (L.GT.NC) THEN
          L = L - NC
          IA(L) = -IA(L)
      ELSE IF (L.GT.0) THEN
          IC(L) = -IC(L)
      ELSE
          L = -L
          IX(L) = -IX(L)
      END IF
      RETURN
      END
* SUBROUTINE PLRMR0               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TRIANGULAR DECOMPOSITION OF KERNEL OF THE ORTHOGONAL PROJECTION IS
* UPDATED AFTER CONSTRAINT DELETION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RA  G(NF)  AUXILIARY VECTOR.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  IOLD  INDEX OF THE OLD ACTIVE CONSTRAINT.
*  IO  KREM  AUXILIARY VARIABLE.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVORT  DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR
*         PLANE ROTATION.
*  S   MXVROT  PLANE ROTATION OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLRMR0(NF,ICA,CR,G,N,IOLD,KREM,IER)
      INTEGER          IER,IOLD,KREM,N,NF
      DOUBLE PRECISION CR(*),G(*)
      INTEGER          ICA(*)
      DOUBLE PRECISION CK,CL
      INTEGER          I,J,K,KC,L,NCA
      NCA = NF - N
      IF (IOLD.LT.NCA) THEN
          K = IOLD* (IOLD-1)/2
          KC = ICA(IOLD)
          CALL MXVCOP(IOLD,CR(K+1),G)
          CALL MXVSET(NCA-IOLD,0.0D0,G(IOLD+1))
          K = K + IOLD
          DO 20 I = IOLD + 1,NCA
              K = K + I
              CALL MXVORT(CR(K-1),CR(K),CK,CL,IER)
              CALL MXVROT(G(I-1),G(I),CK,CL,IER)
              L = K
              DO 10 J = I,NCA - 1
                  L = L + J
                  CALL MXVROT(CR(L-1),CR(L),CK,CL,IER)
   10         CONTINUE
   20     CONTINUE
          K = IOLD* (IOLD-1)/2
          DO 30 I = IOLD,NCA - 1
              L = K + I
              ICA(I) = ICA(I+1)
              CALL MXVCOP(I,CR(L+1),CR(K+1))
              K = L
   30     CONTINUE
          ICA(NCA) = KC
          CALL MXVCOP(NCA,G,CR(K+1))
      END IF
      KREM = 1
      RETURN
      END
* SUBROUTINE PLSETC             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF INITIAL VALUES OF THE CONSTRAINT FUNCTIONS.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CURRENT LINEAR CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  XO(NF)  SAVED VECTOR OF VARIABLES.
*  RU  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CG(NF*MCL)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RA  S(NF)  AUXILIARY VECTOR.
*
* SUBPROGRAMS USED :
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*
      SUBROUTINE PLSETC(NF,NC,X,XO,CF,IC,CG,S)
      INTEGER NF,NC,IC(*)
      DOUBLE PRECISION X(*),XO(*),CF(*),CG(*),S(*)
      DOUBLE PRECISION MXVDOT
      INTEGER JCG,KC
      CALL MXVDIF(NF,X,XO,S)
      JCG=0
      DO 1 KC=1,NC
      IF (IC(KC).NE.0) CF(KC)=CF(KC)+MXVDOT(NF,S,CG(JCG+1))
      JCG=JCG+NF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PLSETG             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* GRADIENT DETERMINATION IN THE FIRST PHASE OF LP SUBROUTINE.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*
* SUBPROGRAMS USED :
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLSETG(NF,NC,IC,CG,G,INEW)
      INTEGER NF,NC,IC(*),INEW
      DOUBLE PRECISION CG(*),G(*)
      INTEGER KC
      CALL MXVSET(NF,0.0D 0,G)
      INEW=0
      DO 4 KC=1,NC
      IF (IC(KC).GE.-10) THEN
      ELSE IF (IC(KC).EQ.-11.OR.IC(KC).EQ.-13.OR.IC(KC).EQ.-15) THEN
      CALL MXVDIR(NF,-1.0D 0,CG((KC-1)*NF+1),G,G)
      INEW=1
      ELSE IF (IC(KC).EQ.-12.OR.IC(KC).EQ.-14.OR.IC(KC).EQ.-16) THEN
      CALL MXVDIR(NF, 1.0D 0,CG((KC-1)*NF+1),G,G)
      INEW=1
      ENDIF
    4 CONTINUE
      RETURN
      END
* SUBROUTINE PLTLAG               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER IS
* COMPUTED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IA(NA)  VECTOR CONTAINING TYPES OF DEVIATIONS.
*  II  IAA(NF+1)  VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS.
*  RI  AZ(NF+1)  VECTOR OF LAGRANGE MULTIPLIERS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  EPS7  TOLERANCE FOR LINEAR AND QUADRATIC PROGRAMMING.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  IO  IOLD  INDEX OF THE REMOVED CONSTRAINT.
*
      SUBROUTINE PLTLAG(NF,N,NC,IX,IA,IAA,AZ,IC,EPS7,UMAX,IOLD)
      INTEGER NF,N,NC,IX(*),IA(*),IAA(*),IC(*),IOLD
      DOUBLE PRECISION AZ(*),EPS7,UMAX
      DOUBLE PRECISION TEMP
      INTEGER NAA,J,K,L
      IOLD=0
      UMAX=0.0D 0
      NAA=NF-N
      DO 2 J=1,NAA
      TEMP=AZ(J)
      L=IAA(J)
      IF (L.GT.NC) THEN
      L=L-NC
      K=IA(L)
      ELSE IF (L.GT.0) THEN
      K=IC(L)
      ELSE
      L=-L
      K=IX(L)
      ENDIF
      IF (K.LE.-5) THEN
      ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN
      ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN
      ELSE
      IOLD=J
      UMAX=ABS(TEMP)
      ENDIF
    2 CONTINUE
      IF (UMAX.LE.EPS7) IOLD=0
      RETURN
      END
* SUBROUTINE PLTRBG               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. LAGRANGE
* MULTIPLIERS ARE DETERMINED. TEST VALUES GMAX AND UMAX ARE COMPUTED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CURRENT LINEAR CONSTRAINTS.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE. VECTOR CZ(1,NF) CONTAINS LAGRANGE
*         MULTIPLIERS BEING DETERMINED.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GN(NF)  TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION IF IT IS
*         NONZERO.
*  RI  EPS7  TOLERANCE FOR LINEAR AND QUADRATIC PROGRAMMING.
*  RO  GMAX  NORM OF THE TRANSFORMED GRADIENT.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  IO  IOLD  INDEX OF THE REMOVED CONSTRAINT.
*
* SUBPROGRAMS USED :
*  S   PLVLAG  GRADIENT IS PREMULTIPLIED BY THE MATRIX WHOSE COLUMNS
*         ARE NORMALS OF THE ACTIVE CONSTRAINTS.
*  S   PLTLAG  COMPUTATION OF THE MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE
*         LAGRANGE MULTIPLIER.
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDPRB  BACK SUBSTITUTION AFTER A CHOLESKI DECOMPOSITION.
*  RF  MXVMAX  L-INFINITY NORM OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,G,GN,EPS7,GMAX,UMAX,
     & IOLD)
      INTEGER NF,N,NC,IX(*),IC(*),ICA(*),IOLD
      DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION MXVMAX
      INTEGER NCA,NCZ
      GMAX=0.0D 0
      IF (N.GT.0) THEN
      CALL MXDRMM(NF,N,CZ,G,GN)
      GMAX=MXVMAX(N,GN)
      ENDIF
      IF (NF.GT.N.AND.GMAX.LE.EPS7) THEN
      NCA=NF-N
      NCZ=N*NF
      CALL PLVLAG(NF,N,NC,ICA,CG,CG,G,CZ(NCZ+1))
      CALL MXDPRB(NCA,CR,CZ(NCZ+1),0)
      CALL PLTLAG(NF,N,NC,IX,IC,ICA,CZ(NCZ+1),IC,EPS7,UMAX,IOLD)
      IF (UMAX.LE.EPS7) IOLD=0
      CALL MXVSET(N,0.0D 0,GN)
      GMAX=0.0D 0
      ELSE
      IOLD=0
      UMAX=0.0D 0
      ENDIF
      RETURN
      END
* SUBROUTINE PLVLAG               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* GRADIENT OF THE OBJECTIVE FUNCTION IS PREMULTIPLIED BY TRANSPOSE
* OF THE MATRIX WHOSE COLUMNS ARE NORMALS OF CURRENT ACTIVE CONSTRAINTS
* AND GRADIENTS OF CURRENT ACTIVE FUNCTIONS.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  II  IAA(NF+1)  VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS.
*  RI  AG(NF*NA)  VECTOR CONTAINING SCALING PARAMETERS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GN(NF+1)  OUTPUT VECTOR.
*
* SUBPROGRAMS USED :
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLVLAG(NF,N,NC,IAA,AG,CG,G,GN)
      INTEGER NF,N,NC,IAA(*)
      DOUBLE PRECISION AG(*),CG(*),G(*),GN(*)
      DOUBLE PRECISION MXVDOT
      INTEGER NAA,J,L
      NAA=NF-N
      DO 1 J=1,NAA
      L=IAA(J)
      IF (L.GT.NC) THEN
      L=L-NC
      GN(J)=MXVDOT(NF,AG((L-1)*NF+1),G)
      ELSE IF (L.GT.0) THEN
      GN(J)=MXVDOT(NF,CG((L-1)*NF+1),G)
      ELSE
      L=-L
      GN(J)=G(L)
      ENDIF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PNINT1                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITH DIRECTIONAL
* DERIVATIVES.
*
* PARAMETERS :
*  RI  RL  LOWER VALUE OF THE STEPSIZE PARAMETER.
*  RI  RU  UPPER VALUE OF THE STEPSIZE PARAMETER.
*  RI  FL  VALUE OF THE OBJECTIVE FUNCTION FOR R=RL.
*  RI  FU  VALUE OF THE OBJECTIVE FUNCTION FOR R=RU.
*  RI  PL  DIRECTIONAL DERIVATIVE FOR R=RL.
*  RI  PU  DIRECTIONAL DERIVATIVE FOR R=RU.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER OBTAINED.
*  II  MODE  MODE OF LINE SEARCH.
*  II  MTYP  METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-QUADRATIC
*         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
*         MTYP=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
*         DERIVATIVES). MTYP=4-CUBIC INTERPOLATION. MTYP=5-CONIC
*         INTERPOLATION.
*  IO  MERR  ERROR INDICATOR. MERR=0 FOR NORMAL RETURN.
*
* METHOD :
* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS.
*
      SUBROUTINE PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR)
      DOUBLE PRECISION RL, RU, FL, FU, PL, PU, R
      INTEGER MODE,MTYP,MERR,NTYP
      DOUBLE PRECISION A,B,C,D,DIS,DEN
      DOUBLE PRECISION C1L,C1U,C2L,C2U,C3L
      PARAMETER (C1L=1.1D 0,C1U=1.0D 3,C2L=1.0D-2,C2U=0.9D 0,
     & C3L=0.1D 0)
      MERR=0
      IF (MODE.LE.0) RETURN
      IF (PL.GE.0.0D 0) THEN
      MERR=2
      RETURN
      ELSE IF (RU.LE.RL) THEN
      MERR=3
      RETURN
      ENDIF
      DO 1 NTYP=MTYP,1,-1
      IF (NTYP.EQ.1) THEN
C
C     BISECTION
C
      IF (MODE.EQ.1) THEN
      R=4.0D 0*RU
      RETURN
      ELSE
      R=0.5D 0*(RL+RU)
      RETURN
      ENDIF
      ELSE IF (NTYP.EQ.MTYP) THEN
      A = (FU-FL)/(PL*(RU-RL))
      B = PU/PL
      ENDIF
      IF (NTYP.EQ.2) THEN
C
C     QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL
C     DERIVATIVE
C
      DEN = 2.0D 0*(1.0D 0-A)
      ELSE IF (NTYP.EQ.3) THEN
C
C     QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL
C     DERIVATIVES
C
      DEN = 1.0D 0 - B
      ELSE IF (NTYP.EQ.4) THEN
C
C     CUBIC EXTRAPOLATION OR INTERPOLATION
C
      C = B - 2.0D 0*A + 1.0D 0
      D = B - 3.0D 0*A + 2.0D 0
      DIS = D*D - 3.0D 0*C
      IF (DIS.LT.0.0D 0) GO TO 1
      DEN = D + SQRT(DIS)
      ELSE IF (NTYP.EQ.5) THEN
C
C     CONIC EXTRAPOLATION OR INTERPOLATION
C
      DIS = A*A - B
      IF (DIS.LT.0.0D 0) GO TO 1
      DEN = A + SQRT(DIS)
      IF (DEN.LE.0.0D 0) GO TO 1
      DEN = 1.0D 0 - B*(1.0D 0/DEN)**3
      ENDIF
      IF (MODE.EQ.1.AND.DEN.GT.0.0D 0.AND.DEN.LT.1.0D 0) THEN
C
C     EXTRAPOLATION ACCEPTED
C
      R = RL + (RU-RL)/DEN
      R = MAX(R,C1L*RU)
      R = MIN(R,C1U*RU)
      RETURN
      ELSE IF (MODE.EQ.2.AND.DEN.GT.1.0D 0) THEN
C
C     INTERPOLATION ACCEPTED
C
      R = RL + (RU-RL)/DEN
      IF (RL.EQ.0.0D 0) THEN
      R = MAX(R,RL+C2L*(RU-RL))
      ELSE
      R = MAX(R,RL+C3L*(RU-RL))
      ENDIF
      R = MIN(R,RL+C2U*(RU-RL))
      RETURN
      ENDIF
    1 CONTINUE
      END
* SUBROUTINE PNINT3                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITHOUT DIRECTIONAL
* DERIVATIVES.
*
* PARAMETERS :
*  RI  RO  INITIAL VALUE OF THE STEPSIZE PARAMETER.
*  RI  RL  LOWER VALUE OF THE STEPSIZE PARAMETER.
*  RI  RU  UPPER VALUE OF THE STEPSIZE PARAMETER.
*  RI  RI  INNER VALUE OF THE STEPSIZE PARAMETER.
*  RI  FO  VALUE OF THE OBJECTIVE FUNCTION FOR R=RO.
*  RI  FL  VALUE OF THE OBJECTIVE FUNCTION FOR R=RL.
*  RI  FU  VALUE OF THE OBJECTIVE FUNCTION FOR R=RU.
*  RI  FI  VALUE OF THE OBJECTIVE FUNCTION FOR R=RI.
*  RO  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER OBTAINED.
*  II  MODE  MODE OF LINE SEARCH.
*  II  MTYP  METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-TWO POINT
*         QUADRATIC INTERPOLATION. MTYP=2-THREE POINT QUADRATIC
*         INTERPOLATION.
*  IO  MERR  ERROR INDICATOR. MERR=0 FOR NORMAL RETURN.
*
* METHOD :
* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS.
*
      SUBROUTINE PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR)
      DOUBLE PRECISION ZERO,HALF,ONE,TWO,THREE,C1L,C1U,C2L,C2U,C3L
      PARAMETER        (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0,
     +                 THREE=3.0D0,C1L=1.1D0,C1U=1.0D3,C2L=1.0D-2,
     +                 C2U=0.9D0,C3L=1.0D-1)
      DOUBLE PRECISION FI,FL,FO,FU,PO,R,RI,RL,RO,RU
      INTEGER          MERR,MODE,MTYP
      DOUBLE PRECISION AI,AL,AU,DEN,DIS
      INTEGER          NTYP
      LOGICAL          L1,L2
      MERR = 0
      IF (MODE.LE.0) RETURN
      IF (PO.GE.ZERO) THEN
          MERR = 2
          RETURN

      ELSE IF (RU.LE.RL) THEN
          MERR = 3
          RETURN
      END IF
      L1 = RL .LE. RO
      L2 = RI .LE. RL
      DO 10 NTYP = MTYP,1,-1
          IF (NTYP.EQ.1) THEN
*
*     BISECTION
*
              IF (MODE.EQ.1) THEN
                  R = TWO*RU
                  RETURN
              ELSE IF (RI-RL.LE.RU-RI) THEN
                  R = HALF* (RI+RU)
                  RETURN
              ELSE
                  R = HALF* (RL+RI)
                  RETURN
              END IF
          ELSE IF (NTYP.EQ.MTYP .AND. L1) THEN
              IF (.NOT.L2) AI = (FI-FO)/ (RI*PO)
              AU = (FU-FO)/ (RU*PO)
          END IF
          IF (L1 .AND. (NTYP.EQ.2.OR.L2)) THEN
*
*     TWO POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION
*
              IF (AU.GE.ONE) GO TO 10
              R = HALF*RU/ (ONE-AU)
          ELSE IF (.NOT.L1 .OR. .NOT.L2 .AND. NTYP.EQ.3) THEN
*
*     THREE POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION
*
              AL = (FI-FL)/ (RI-RL)
              AU = (FU-FI)/ (RU-RI)
              DEN = AU - AL
              IF (DEN.LE.ZERO) GO TO 10
              R = RI - HALF* (AU* (RI-RL)+AL* (RU-RI))/DEN
          ELSE IF (L1 .AND. .NOT.L2 .AND. NTYP.EQ.4) THEN
*
*     THREE POINT CUBIC EXTRAPOLATION OR INTERPOLATION
*
              DIS = (AI-ONE)* (RU/RI)
              DEN = (AU-ONE)* (RI/RU) - DIS
              DIS = AU + AI - DEN - TWO* (ONE+DIS)
              DIS = DEN*DEN - THREE*DIS
              IF (DIS.LT.ZERO) GO TO 10
              DEN = DEN + SQRT(DIS)
              IF (DEN.EQ.ZERO) GO TO 10
              R = (RU-RI)/DEN
          ELSE
              GO TO 10
          END IF
          IF (MODE.EQ.1 .AND. R.GT.RU) THEN
*
*     EXTRAPOLATION ACCEPTED
*
              R = MAX(R,C1L*RU)
              R = MIN(R,C1U*RU)
              RETURN
          ELSE IF (MODE.EQ.2 .AND. R.GT.RL .AND. R.LT.RU) THEN
*
*     INTERPOLATION ACCEPTED
*
              IF (RI.EQ.ZERO .AND. NTYP.NE.4) THEN
                  R = MAX(R,RL+C2L* (RU-RL))
              ELSE
                  R = MAX(R,RL+C3L* (RU-RL))
              END IF
              R = MIN(R,RL+C2U* (RU-RL))
              IF (R.EQ.RI) GO TO 10
              RETURN
          END IF
   10 CONTINUE
      END
* SUBROUTINE PNSTEP                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/01/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF A SCALING FACTOR FOR THE BOUNDARY STEP.
*
* PARAMETERS :
*  RI  DEL  MAXIMUM STEPSIZE.
*  RI  A  INPUT PARAMETER.
*  RI  B  INPUT PARAMETER.
*  RI  C  INPUT PARAMETER.
*  RO  ALF  SCALING FACTOR FOR THE BOUNDARY STEP SUCH THAT
*         A**2+2*B*ALF+C*ALF**2=DEL**2.
*
      SUBROUTINE PNSTEP(DEL,A,B,C,ALF)
      DOUBLE PRECISION  DEL, A, B, C, ALF
      DOUBLE PRECISION  DEN, DIS
      DOUBLE PRECISION  ZERO
      PARAMETER (ZERO = 0.0D 0)
      ALF = ZERO
      DEN = (DEL+A) * (DEL-A)
      IF (DEN .LE. ZERO) RETURN
      DIS = B*B + C*DEN
      IF (B .GE. ZERO) THEN
      ALF = DEN / (SQRT(DIS) + B)
      ELSE
      ALF = (SQRT(DIS) - B) / C
      ENDIF
      RETURN
      END
* SUBROUTINE PP0AF8             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF VALUE OF THE AUGMENTED LAGRANGIAN FUNCTION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  N  DIMENSION OF THE CONSTRAINT NULL SPACE.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  CF(NC+1)  VECTOR CONTAINING VALUES OF THE CONSTRAINTS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CZ(NC)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  RPF  PENALTY COEFFICIENT.
*  RO  FC  VALUE OF THE PENALTY TERM.
*  RO  F  VALUE OF THE PENALTY FUNCTION.
*
      SUBROUTINE PP0AF8(NF,N,NC,CF,IC,ICA,CL,CU,CZ,RPF,FC,F)
      INTEGER NF,N,NC,IC(*),ICA(*)
      DOUBLE PRECISION CF(*),CL(*),CU(*),CZ(*),RPF,FC,F
      DOUBLE PRECISION POM,TEMP
      INTEGER J,KC
      FC=0.0D0
      DO 1 KC=1,NC
      IF (IC(KC).GT.0) THEN
      POM=0.0D0
      TEMP=CF(KC)
      IF (IC(KC).EQ.1.OR.IC(KC).GE.3) POM=MIN(POM,TEMP-CL(KC))
      IF (IC(KC).EQ.2.OR.IC(KC).GE.3) POM=MIN(POM,CU(KC)-TEMP)
      FC=FC+RPF*ABS(POM)
      ENDIF
    1 CONTINUE
      DO 2 J=1,NF-N
      KC=ICA(J)
      IF (KC.GT.0) THEN
      POM=0.0D0
      TEMP=CF(KC)
      IF (IC(KC).EQ.1.OR.IC(KC).EQ.3.OR.IC(KC).EQ.5)
     & POM=MIN(POM,TEMP-CL(KC))
      IF (IC(KC).EQ.2.OR.IC(KC).EQ.4.OR.IC(KC).EQ.6)
     & POM=MAX(POM,TEMP-CU(KC))
      FC=FC-CZ(J)*POM
      ENDIF
    2 CONTINUE
      F=CF(NC+1)+FC
      RETURN
      END
* SUBROUTINE PPSET2             ALL SYSTEMS                   97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE NEW PENALTY PARAMETERS.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CZ(NF)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  CP(NC)  VECTOR CONTAINING PENALTY PARAMETERS.
*
      SUBROUTINE PPSET2(NF,N,NC,ICA,CZ,CP)
      INTEGER NF,N,NC,ICA(*)
      DOUBLE PRECISION CZ(*),CP(*)
      DOUBLE PRECISION TEMP
      INTEGER J,L,KC
      DO 1 KC=1,NC
      CP(KC)=0.5D 0*CP(KC)
    1 CONTINUE
      DO 2 J=1,NF-N
      L=ICA(J)
      IF (L.GT.0) THEN
      TEMP=ABS(CZ(J))
      CP(L)=MAX(TEMP,CP(L)+0.5D 0*TEMP)
      ENDIF
    2 CONTINUE
      RETURN
      END
* SUBROUTINE PS0G01                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SIMPLE SEARCH WITH TRUST REGION UPDATE.
*
* PARAMETERS :
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PP  QUADRATIC PART OF THE PREDICTED FUNCTION VALUE.
*  RU  XDEL  TRUST REGION BOUND.
*  RO  XDELO  PREVIOUS TRUST REGION BOUND.
*  RI  XMAX MAXIMUM STEPSIZE.
*  RI  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
*  RI  SNORM  EUCLIDEAN NORM OF THE DIRECTION VECTOR.
*  RI  BET1  LOWER BOUND FOR STEPSIZE REDUCTION.
*  RI  BET2  UPPER BOUND FOR STEPSIZE REDUCTION.
*  RI  GAM1  LOWER BOUND FOR STEPSIZE EXPANSION.
*  RI  GAM2  UPPER BOUND FOR STEPSIZE EXPANSION.
*  RI  EPS4  FIRST TOLERANCE FOR RATIO DF/DFPRED. STEP BOUND IS
*         DECREASED IF DF/DFPRED<EPS4.
*  RI  EPS5  SECOND TOLERANCE FOR RATIO DF/DFPRED. STEP BOUND IS
*         INCREASED IF IT IS ACTIVE AND DF/DFPRED>EPS5.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  IU  IDIR INDICATOR FOR DIRECTION DETERMINATION.
*         IDIR=0-BASIC DETERMINATION. IDIR=1-DETERMINATION
*         AFTER STEPSIZE REDUCTION. IDIR=2-DETERMINATION AFTER
*         STEPSIZE EXPANSION.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-STEP
*         BOUND WAS DECREASED. ITERS=2-STEP BOUND WAS UNCHANGED.
*         ITERS=3-STEP BOUND WAS INCREASED. ITERS=6-FIRST STEPSIZE.
*  II  ITERD TERMINATION INDICATOR. ITERD<0-BAD DECOMPOSITION.
*         ITERD=0-DESCENT DIRECTION. ITERD=1-NEWTON LIKE STEP.
*         ITERD=2-INEXACT NEWTON LIKE STEP. ITERD=3-BOUNDARY STEP.
*         ITERD=4-DIRECTION WITH THE NEGATIVE CURVATURE.
*         ITERD=5-MARQUARDT STEP.
*  IO  MAXST  MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
*         STEPSIZE WAS NOT OR WAS REACHED.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  MRED  MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  KTERS  TERMINATION SELECTION. KTERS=1-NORMAL TERMINATION.
*         KTERS=6-FIRST STEPSIZE.
*  II  MES1  SWITCH FOR EXTRAPOLATION. MES1=1-CONSTANT INCREASING OF
*         THE INTERVAL. MES1=2-EXTRAPOLATION SPECIFIED BY THE PARAMETER
*         MES. MES1=3 SUPPRESSED EXTRAPOLATION.
*  II  MES2  SWITCH FOR TERMINATION. MES2=1-NORMAL TERMINATION.
*         MES2=2-TERMINATION AFTER AT LEAST TWO STEPS (ASYMPTOTICALLY
*         PERFECT LINE SEARCH).
*  II  MES3  SAFEGUARD AGAINST ROUNDING ERRORS. MES3=0-SAFEGUARD
*         SUPPRESSED. MES3=1-FIRST LEVEL OF SAFEGUARD. MES3=2-SECOND
*         LEVEL OF SAFEGUARD.
*  IU  ISYS  CONTROL PARAMETER.
*
* COMMON DATA :
*
* METHOD :
* G.A.SCHULTZ, R.B.SCHNABEL, R.H.BYRD: A FAMILY OF TRUST-REGION-BASED
* ALGORITHMS FOR UNCONSTRAINED MINIMIZATION WITH STRONG GLOBAL
* CONVERGENCE PROPERTIES, SIAM J. NUMER.ANAL. 22 (1985) PP. 47-67.
*
      SUBROUTINE PS0G01(R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,
     & BET1,BET2,GAM1,GAM2,EPS4,EPS5,KD,LD,IDIR,ITERS,ITERD,
     & MAXST,NRED,MRED,KTERS,MES1,MES2,MES3,ISYS)
      INTEGER KD,LD,IDIR,ITERS,ITERD,MAXST,NRED,MRED,KTERS,
     & MES1,MES2,MES3,ISYS
      DOUBLE PRECISION R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,BET1,
     & BET2,GAM1,GAM2,EPS4,EPS5
      DOUBLE PRECISION DF,DFPRED
      INTEGER NRED1,NRED2
      SAVE NRED1,NRED2
      IF (ISYS.EQ.1) GO TO 2
C      GO TO (1,2) ISYS+1
    1 CONTINUE
      IF (IDIR.EQ.0) THEN
      NRED1=0
      NRED2=0
      ENDIF
      IDIR=0
      XDELO=XDEL
C
C     COMPUTATION OF THE NEW FUNCTION VALUE
C
      R=MIN(1.0D 0,RMAX)
      KD= 0
      LD=-1
      ISYS=1
      RETURN
    2 CONTINUE
      IF(KTERS.LT.0.OR.KTERS.GT.5) THEN
      ITERS=6
      ELSE
      DF=FO-F
      DFPRED=-R*(PO+R*PP)
      IF (DF.LT.EPS4*DFPRED) THEN
C
C     STEP IS TOO LARGE, IT HAS TO BE REDUCED
C
      IF (MES1.EQ.1) THEN
      XDEL=BET2*SNORM
      ELSE IF (MES1.EQ.2) THEN
      XDEL=BET2*MIN(0.5D 0*XDEL,SNORM)
      ELSE
      XDEL=0.5D 0*PO*SNORM/(PO+DF)
      XDEL=MAX(XDEL,BET1*SNORM)
      XDEL=MIN(XDEL,BET2*SNORM)
      ENDIF
      ITERS=1
      IF (MES3.LE.1) THEN
      NRED2=NRED2+1
      ELSE
      IF (ITERD.GT.2) NRED2=NRED2+1
      ENDIF
      ELSE IF (DF.LE.EPS5*DFPRED) THEN
C
C     STEP IS SUITABLE
C
      ITERS=2
      ELSE
C
C     STEP IS TOO SMALL, IT HAS TO BE ENLARGED
C
      IF (MES2.EQ.2) THEN
      XDEL=MAX(XDEL,GAM1*SNORM)
      ELSE IF (ITERD.GT.2) THEN
      XDEL=GAM1*XDEL
      ENDIF
      ITERS=3
      ENDIF
      XDEL=MIN(XDEL,XMAX,GAM2*SNORM)
      IF (FO.LE.F) THEN
      IF (NRED1.GE.MRED) THEN
      ITERS=-1
      ELSE
      IDIR=1
      ITERS=0
      NRED1=NRED1+1
      ENDIF
      ENDIF
      ENDIF
      MAXST=0
      IF (XDEL.GE.XMAX) MAXST=1
      IF (MES3.EQ.0) THEN
      NRED=NRED1
      ELSE
      NRED=NRED2
      ENDIF
      ISYS=0
      RETURN
      END
* SUBROUTINE PS0L02                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
*  EXTENDED LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES.
*
* PARAMETERS :
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  RO  INITIAL VALUE OF THE STEPSIZE PARAMETER.
*  RO  RP  PREVIOUS VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  PP  PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  FMIN  LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FMAX  UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  RMIN  MINIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  TOLS  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
*         CHANGE OF THE FUNCTION VALUE).
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  MRED  MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  IO  MAXST  MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
*         STEPSIZE WAS NOT OR WAS REACHED.
*  II  IEST  LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND
*         IS NOT OR IS GIVEN.
*  II  INITS  CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE
*         IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL
*         STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED
*         INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL
*         STEPSIZE.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
*         LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
*         STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
*         ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
*         ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
*         ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
*         DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
*  II  KTERS  TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH.
*         KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE.
*         KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE.
*         KTERS=6-FIRST STEPSIZE.
*  II  MES  METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC
*         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
*         MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
*         DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC
*         INTERPOLATION.
*  IU  ISYS  CONTROL PARAMETER.
*
* SUBPROGRAM USED :
*  S   PNINT3  EXTRAPOLATION OR INTERPOLATION WITHOUT DIRECTIONAL
*         DERIVATIVES.
*
* METHOD :
* SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH EXTENDED TERMINATION
* CRITERIA.
*
      SUBROUTINE PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,
     & TOLS,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS)
      INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS
      DOUBLE PRECISION R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS
      DOUBLE PRECISION RL,FL,RU,FU,RI,FI,RTEMP,TOL
      INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2
      LOGICAL L1,L2,L3,L4,L6,L7
      PARAMETER(TOL=1.0D-4)
      SAVE MTYP,MODE,MES1,MES2
      SAVE RL,FL,RU,FU,RI,FI
      IF (ISYS.EQ.1) GO TO 3
C      GO TO (1,3) ISYS+1
    1 CONTINUE
      MES1=2
      MES2=2
      ITERS=0
      IF (PO.GE.0.0D 0) THEN
      R=0.0D 0
      ITERS=-2
      GO TO 4
      ENDIF
      IF(RMAX.LE.0.0D 0) THEN
      ITERS= 0
      GO TO 4
      ENDIF
C
C     INITIAL STEPSIZE SELECTION
C
      IF (INITS.GT.0) THEN
      RTEMP=FMIN-F
      ELSE IF (IEST.EQ.0) THEN
      RTEMP=F-FP
      ELSE
      RTEMP=MAX(F-FP,1.0D 1*(FMIN-F))
      ENDIF
      INIT1=ABS(INITS)
      RP=0.0D 0
      FP=FO
      PP=PO
      IF (INIT1.EQ.0) THEN
      ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN
      R=1.0D 0
      ELSE IF (INIT1.EQ.2) THEN
      R=MIN(1.0D 0,4.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.3) THEN
      R=MIN(1.0D 0, 2.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.4) THEN
      R=2.0D 0*RTEMP/PO
      ENDIF
      RTEMP=R
      R=MAX(R,RMIN)
      R=MIN(R,RMAX)
      MODE=0
      RL=0.0D 0
      FL=FO
      RU=0.0D 0
      FU=FO
      RI=0.0D 0
      FI=FO
C
C     NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION)
C
    2 CALL PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR)
      IF (MERR.GT.0) THEN
      ITERS=-MERR
      GO TO 4
      ELSE IF (MODE.EQ.1) THEN
      NRED=NRED-1
      R=MIN(R,RMAX)
      ELSE IF (MODE.EQ.2) THEN
      NRED=NRED+1
      ENDIF
C
C     COMPUTATION OF THE NEW FUNCTION VALUE
C
      KD= 0
      LD=-1
      ISYS=1
      RETURN
    3 CONTINUE
      IF (ITERS.NE.0) GO TO 4
      IF (F.LE.FMIN) THEN
      ITERS=7
      GO TO 4
      ELSE
      L1=R.LE.RMIN.AND.NIT.NE.KIT
      L2=R.GE.RMAX
      L3=F-FO.LE.TOLS*R*PO.OR.F-FMIN.LE.(FO-FMIN)/1.0D 1
      L4=F-FO.GE.(1.0D 0-TOLS)*R*PO.OR.MES2.EQ.2.AND.MODE.EQ.2
      L6=RU-RL.LE.TOL*RU.AND.MODE.EQ.2
      L7=MES2.LE.2.OR.MODE.NE.0
      MAXST=0
      IF (L2) MAXST=1
      ENDIF
C
C     TEST ON TERMINATION
C
      IF (L1.AND..NOT.L3) THEN
      ITERS=0
      GO TO 4
      ELSE IF (L2.AND..NOT.F.GE.FU) THEN
      ITERS=7
      GO TO 4
      ELSE IF (L6) THEN
      ITERS=1
      GO TO 4
      ELSE IF (L3.AND.L7.AND.KTERS.EQ.5) THEN
      ITERS=5
      GO TO 4
      ELSE IF (L3.AND.L4.AND.L7.AND.(KTERS.EQ.2.OR.KTERS.EQ.3.OR.
     *         KTERS.EQ.4)) THEN
      ITERS=2
      GO TO 4
      ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN
      ITERS=6
      GO TO 4
      ELSE IF(ABS(NRED).GE.MRED) THEN
      ITERS=-1
      GO TO 4
      ELSE
      RP=R
      FP=F
      MODE=MAX(MODE,1)
      MTYP=ABS(MES)
      IF(F.GE.FMAX) MTYP=1
      ENDIF
      IF (MODE.EQ.1) THEN
C
C     INTERVAL CHANGE AFTER EXTRAPOLATION
C
      RL=RI
      FL=FI
      RI=RU
      FI=FU
      RU=R
      FU=F
      IF (F.GE.FI) THEN
      NRED=0
      MODE=2
      ELSE IF ( MES1 .EQ. 1) THEN
      MTYP=1
      ENDIF
C
C     INTERVAL CHANGE AFTER INTERPOLATION
C
      ELSE IF (R.LE.RI) THEN
      IF (F.LE.FI) THEN
      RU=RI
      FU=FI
      RI=R
      FI=F
      ELSE
      RL=R
      FL=F
      ENDIF
      ELSE
      IF (F.LE.FI) THEN
      RL=RI
      FL=FI
      RI=R
      FI=F
      ELSE
      RU=R
      FU=F
      ENDIF
      ENDIF
      GO TO 2
    4 ISYS=0
      RETURN
      END
* SUBROUTINE PS1L01                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
*  STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES.
*
* PARAMETERS :
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  RP  PREVIOUS VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  PP  PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  FMIN  LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FMAX  UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  RMIN  MINIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  TOLS  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
*         CHANGE OF THE FUNCTION VALUE).
*  RI  TOLP  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
*         CHANGE OF THE DIRECTIONAL DERIVATIVE).
*  RO  PAR1  PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
*         UPDATES.
*  RO  PAR2  PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
*         UPDATES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  MRED  MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  IO  MAXST  MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
*         STEPSIZE WAS NOT OR WAS REACHED.
*  II  IEST  LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND
*         IS NOT OR IS GIVEN.
*  II  INITS  CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE
*         IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL
*         STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED
*         INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL
*         STEPSIZE.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
*         LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
*         STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
*         ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
*         ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
*         ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
*         DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
*  II  KTERS  TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH.
*         KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE.
*         KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE.
*         KTERS=6-FIRST STEPSIZE.
*  II  MES  METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC
*         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
*         MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
*         DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC
*         INTERPOLATION.
*  IU  ISYS  CONTROL PARAMETER.
*
* SUBPROGRAM USED :
*  S   PNINT1  EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL
*         DERIVATIVES.
*
* METHOD :
* SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION
* CRITERIA.
*
      SUBROUTINE PS1L01(R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX,
     & TOLS,TOLP,PAR1,PAR2,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,
     & ITERS,KTERS,MES,ISYS)
      INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS
      DOUBLE PRECISION R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX,
     & TOLS,TOLP,PAR1,PAR2
      DOUBLE PRECISION RL,FL,PL,RU,FU,PU,RTEMP
      INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2,MES3
      LOGICAL L1,L2,L3,L5,L7,M1,M2,M3
      DOUBLE PRECISION CON,CON1
      PARAMETER (CON=1.0D-2,CON1=1.0D-13)
      SAVE MTYP,MODE,MES1,MES2,MES3
      SAVE RL,FL,PL,RU,FU,PU
      IF (ISYS.EQ.1) GO TO 3
C      GO TO (1,3) ISYS+1
    1 CONTINUE
      MES1=2
      MES2=2
      MES3=2
      ITERS=0
      IF (PO.GE.0.0D 0) THEN
      R=0.0D 0
      ITERS=-2
      GO TO 4
      ENDIF
      IF(RMAX.LE.0.0D 0) THEN
      ITERS=0
      GO TO 4
      ENDIF
C
C     INITIAL STEPSIZE SELECTION
C
      IF (INITS.GT.0) THEN
      RTEMP=FMIN-F
      ELSE IF (IEST.EQ.0) THEN
      RTEMP=F-FP
      ELSE
      RTEMP=MAX(F-FP,1.0D 1*(FMIN-F))
      ENDIF
      INIT1=ABS(INITS)
      RP=0.0D 0
      FP=FO
      PP=PO
      IF (INIT1.EQ.0) THEN
      ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN
      R=1.0D 0
      ELSE IF (INIT1.EQ.2) THEN
      R=MIN(1.0D 0,4.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.3) THEN
      R=MIN(1.0D 0, 2.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.4) THEN
      R=2.0D 0*RTEMP/PO
      ENDIF
      R=MAX(R,RMIN)
      R=MIN(R,RMAX)
      MODE=0
      RU=0.0D 0
      FU=FO
      PU=PO
C
C     NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION)
C
    2 CALL PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR)
      IF (MERR.GT.0) THEN
      ITERS=-MERR
      GO TO 4
      ELSE IF (MODE.EQ.1) THEN
      NRED=NRED-1
      R=MIN(R,RMAX)
      ELSE IF (MODE.EQ.2) THEN
      NRED=NRED+1
      ENDIF
C
C     COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL
C     DERIVATIVE
C
      KD= 1
      LD=-1
      ISYS=1
      RETURN
    3 CONTINUE
      IF (MODE.EQ.0) THEN
      PAR1=P/PO
      PAR2=F-FO
      ENDIF
      IF (ITERS.NE.0) GO TO 4
      IF (F.LE.FMIN) THEN
      ITERS=7
      GO TO 4
      ELSE
      L1=R.LE.RMIN.AND.NIT.NE.KIT
      L2=R.GE.RMAX
      L3=F-FO.LE.TOLS*R*PO
      L5=P.GE.TOLP*PO.OR.MES2.EQ.2.AND.MODE.EQ.2
      L7=MES2.LE.2.OR.MODE.NE.0
      M1=.FALSE.
      M2=.FALSE.
      M3=L3
      IF(MES3.GE.1) THEN
      M1=ABS(P).LE.CON*ABS(PO).AND.FO-F.GE.(CON1/CON)*ABS(FO)
      L3=L3.OR.M1
      ENDIF
      IF(MES3.GE.2) THEN
      M2=ABS(P).LE.0.5D 0*ABS(PO).AND.ABS(FO-F).LE.2.0D 0*CON1*ABS(FO)
      L3=L3.OR.M2
      ENDIF
      MAXST=0
      IF (L2) MAXST=1
      ENDIF
C
C     TEST ON TERMINATION
C
      IF (L1.AND..NOT.L3) THEN
      ITERS=0
      GO TO 4
      ELSE IF (L2.AND.L3.AND..NOT.L5)  THEN
      ITERS=7
      GO TO 4
      ELSE IF (M3.AND.MES1.EQ.3) THEN
      ITERS=5
      GO TO 4
      ELSE IF (L3.AND.L5.AND.L7) THEN
      ITERS=4
      GO TO 4
      ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN
      ITERS=6
      GO TO 4
      ELSE IF (ABS(NRED).GE.MRED) THEN
      ITERS=-1
      GO TO 4
      ELSE
      RP=R
      FP=F
      PP=P
      MODE=MAX(MODE,1)
      MTYP=ABS(MES)
      IF(F.GE.FMAX) MTYP=1
      ENDIF
      IF (MODE.EQ.1) THEN
C
C     INTERVAL CHANGE AFTER EXTRAPOLATION
C
      RL=RU
      FL=FU
      PL=PU
      RU=R
      FU=F
      PU=P
      IF (.NOT.L3) THEN
      NRED=0
      MODE=2
      ELSE IF ( MES1 .EQ. 1) THEN
      MTYP=1
      ENDIF
      ELSE
C
C     INTERVAL CHANGE AFTER INTERPOLATION
C
      IF (.NOT.L3) THEN
      RU=R
      FU=F
      PU=P
      ELSE
      RL=R
      FL=F
      PL=P
      ENDIF
      ENDIF
      GO TO 2
    4 ISYS=0
      RETURN
      END
* SUBROUTINE PUDBG1                ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
* USING THE FACTORIZATION B=L*D*TRANS(L).
*
* PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RU  H(M)  FACTORIZATION B=L*D*TRANS(L) OF A POSITIVE
*         DEFINITE APPROXIMATION OF THE HESSIAN MATRIX.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  GO(NF)  GRADIENTS DIFFERENCE.
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  IO  ITERH  TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
*         ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
*  II  MET1  SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED.
*         MET1=2 INITIAL SELF SCALING.
*  II  MEC  CORRECTION IF THE NEGATIVE CURVATURE OCCURS.
*         MEC=1-CORRECTION SUPPRESSED. MEC=2-POWELL'S CORRECTION.
*
* SUBPROGRAMS USED :
*  S   MXDPGU  CORRECTION OF A DENSE SYMMETRIC POSITIVE DEFINITE
*         MATRIX IN THE FACTORED FORM B=L*D*TRANS(L).
*  S   MXDPGS  SCALING OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
*         IN THE FACTORED FORM B=L*D*TRANS(L).
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  RF  MXVDOT  DOT PRODUCT OF VECTORS.
*  S   MXVSCL  SCALING OF A VECTOR.
*
* METHOD :
* BFGS VARIABLE METRIC METHOD.
*
      SUBROUTINE PUDBG1(N,H,G,S,XO,GO,R,PO,NIT,KIT,ITERH,MET,MET1,MEC)
      DOUBLE PRECISION PO,R
      INTEGER          ITERH,KIT,MET,MET1,MEC,N,NIT
      DOUBLE PRECISION G(*),GO(*),H(*),S(*),XO(*)
      DOUBLE PRECISION A,B,C,GAM,PAR,DEN,DIS
      LOGICAL          L1,L3
      DOUBLE PRECISION MXVDOT,MXDPGP
      L1 = MET1 .GE. 3 .OR. MET1 .EQ. 2 .AND. NIT .EQ. KIT
      L3 = .NOT. L1
*
*     DETERMINATION OF THE PARAMETERS B, C
*
      B = MXVDOT(N,XO,GO)
      A = 0.0D0
      IF (L1) THEN
          CALL MXVCOP(N,GO,S)
          CALL MXDPGB(N,H,S,1)
          A=MXDPGP(N,H,S,S)
          IF (A.LE.0.0D0) THEN
              ITERH=1
              RETURN
          END IF
      END IF
      CALL MXVDIF(N,GO,G,S)
      CALL MXVSCL(N,R,S,S)
      C = -R*PO
      IF (C.LE.0.0D0) THEN
          ITERH = 3
          RETURN
      END IF
      IF (MEC.GT.1) THEN
          IF (B.LE.1.0D-4*C) THEN
*
*     POWELL'S CORRECTION
*
              DIS=(1.0D0-0.1D0)*C/(C-B)
              CALL MXVDIF(N,GO,S,GO)
              CALL MXVDIR(N,DIS,GO,S,GO)
              B=C+DIS*(B-C)
              IF (L1) A=C+2.0D0*(1.0D0-DIS)*(B-C)+DIS*DIS*(A-C)
          ENDIF
      ELSE
          IF (B.LE.1.0D-4*C) THEN
              ITERH = 2
              RETURN
          ENDIF
      ENDIF
      IF (L1) THEN
*
*     DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
*
          IF (MET.EQ.1) THEN
              PAR = C/B
          ELSE IF (A.LE.0.0D 0) THEN
              PAR = C/B
          ELSE
              PAR=SQRT(C/A)
          ENDIF
          GAM = PAR
          IF (MET1.GT.1) THEN
              IF (NIT.NE.KIT) THEN
                  L3=GAM.LT.0.5D0.OR.GAM.GT.4.0D0
              ENDIF
          ENDIF
      ENDIF
      IF (L3) THEN
          GAM = 1.0D0
          PAR = GAM
      END IF
      IF (MET.EQ.1) THEN
*
*     BFGS UPDATE
*
      CALL MXDPGU(N,H,PAR/B,GO,XO)
      CALL MXDPGU(N,H,-1.0D0/C,S,XO)
      ELSE
*
*     HOSHINO UPDATE
*
      DEN=PAR*B+C
      DIS=0.5D0*B
      CALL MXVDIR(N,PAR,GO,S,S)
      CALL MXDPGU(N,H,PAR/DIS,GO,XO)
      CALL MXDPGU(N,H,-1.0D0/DEN,S,XO)
      ENDIF
      ITERH = 0
      IF (GAM.EQ.1.0D0) RETURN
      CALL MXDPGS(N,H,1.0D0/GAM)
      RETURN
      END
* SUBROUTINE PUDBI1                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VARIABLE METRIC UPDATES.
*
* PARAMETERS :
*  II  N  NUMBER OF VARIABLES.
*  RU  H(N*(N+1)/2)  UPDATED APPROXIMATION OF THE INVERSE HESSIAN
*         MATRIX.
*  RA  S(N)  AUXILIARY VECTOR.
*  RI  XO(N)  VECTOR OF VARIABLES DIFFERENCE.
*  RI  GO(N)  GRADIENT DIFFERENCE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  PAR1  PARAMETER FOR CONTROL SCALING.
*  RO  PAR2  PARAMETER FOR CONTROL SCALING.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RI  P  CURRENT VALUE OF THE DIRECTIONAL DERIVATIVE.
*  II  NIT  NUMBER OF ITERATIONS.
*  II  KIT  INDEX OF THE ITERATION WITH THE LAST RESTART.
*  II  MET  VARIABLE METRIC UPDATE.
*  II  MET1  SCALING STRATEGY.
*  II  MET2  CORRECTION RULE.
*  IU  IDECF  DECOMPOSITION INDICATOR.
*  II  ITERD  TERMINATION INDICATOR. ITERD<0-BAD DECOMPOSITION.
*         ITERD=0-DESCENT DIRECTION. ITERD=1-NEWTON LIKE STEP.
*         ITERD=2-INEXACT NEWTON LIKE STEP. ITERD=3-BOUNDARY STEP.
*         ITERD=4-DIRECTION WITH THE NEGATIVE CURVATURE.
*         ITERD=5-MARQUARDT STEP.
*  IO  ITERH  UPDATE INDICATOR. ITERH=0-SUCCESSFUL UPDATE.
*         ITERH>0-UNSUCCESSFUL UPDATE.
*
* METHOD :
* VARIOUS VARIABLE METRIC UPDATES INCLUDING BFGS UPDATE.
*
      SUBROUTINE PUDBI1(N,H,S,XO,GO,R,PO,PAR1,PAR2,F,FO,P,NIT,KIT,
     & MET,MET1,MET2,IDECF,ITERD,ITERH)
      INTEGER N,NIT,KIT,MET,MET1,MET2,IDECF,ITERD,ITERH
      DOUBLE PRECISION H(*),S(*),XO(*),GO(*),R,PO
      DOUBLE PRECISION PAR1,PAR2
      DOUBLE PRECISION F,FO,P
      DOUBLE PRECISION AA,CC
      DOUBLE PRECISION MXVDOT
      DOUBLE PRECISION DIS,POM,POM3,POM4,A,B,C,GAM,RHO,PAR
      DOUBLE PRECISION DEN
      LOGICAL L1,L2,L3
      IF (MET.LE.0) GO TO 22
      IF (IDECF.NE.9) THEN
      ITERH=-1
      GO TO 22
      ENDIF
      L1=ABS(MET1).GE.3.OR.ABS(MET1).EQ.2.AND.NIT.EQ.KIT
      L3=.NOT.L1
*
*     DETERMINATION OF THE PARAMETERS A, B, C
*
      B=MXVDOT(N,XO,GO)
      IF (B.LE.0.0D 0) THEN
      ITERH=2
      GO TO 22
      ENDIF
      CALL MXDSMM(N,H,GO,S)
      A=MXVDOT(N,GO,S)
      IF (A.LE.0.0D 0) THEN
      ITERH=1
      GO TO 22
      ENDIF
      IF(MET.NE.1.OR.L1) THEN
      IF (ITERD.NE.1) THEN
      C=0.0D 0
      ELSE
      C=-R*PO
      IF (C.LE.0.0D 0) THEN
      ITERH=3
      GO TO 22
      ENDIF
      ENDIF
      ELSE
      C=0.0D 0
      ENDIF
*
*     DETERMINATION OF THE PARAMETER RHO (NONQUADRATIC PROPERTIES)
*
      IF (MET2.EQ.1) THEN
      RHO=1.0D 0
      ELSE IF (FO-F+P.EQ.0) THEN
      RHO=1.0D 0
      ELSE
      RHO=0.5D 0*B/(FO-F+P)
      ENDIF
      IF(RHO.LE.1.0D-2) RHO=1.0D 0
      IF(RHO.GE.1.0D 2) RHO=1.0D 0
      AA=A/B
      CC=C/B
      IF (L1) THEN
*
*     DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
*
      PAR=A/B
      POM3=0.7D 0
      POM4=6.0D 0
      GAM=RHO/PAR
      IF (NIT.NE.KIT) THEN
      IF (MET1.EQ.3) THEN
      L2=PAR2.LE.0.0D 0
      L3=L2.AND.ABS(PAR1).LE.0.2D 0
      L3=L3.OR.(.NOT.L2.AND.GAM.GT.1.0D 0)
      L3=L3.OR.(L2.AND.PAR1.LT.0.0D 0.AND.GAM.GT.1.0D 0)
      L3=L3.OR.(L2.AND.PAR1.GT.0.0D 0.AND.GAM.LT.1.0D 0)
      L3=L3.OR.GAM.LT.POM3
      L3=L3.OR.GAM.GT.POM4
      ELSE IF (MET1.EQ.4) THEN
      L3=GAM.LT.POM3.OR.GAM.GT.POM4
      ENDIF
      ENDIF
      ENDIF
      IF (L3) THEN
      GAM=1.0D 0
      PAR=RHO/GAM
      ENDIF
      IF (MET.EQ.1) GO TO 17
*
*     NEW UPDATE
*
      POM=1.0D 0/(AA*CC)
      IF (POM.LT.1.0D 0) THEN
      POM=MAX(1.0D-15,(SQRT(C/A)-POM)/(1.0D 0-POM))
      GO TO 20
      ENDIF
   17 CONTINUE
*
*     BFGS UPDATE
*
      POM=1.0D 0
      DIS=PAR+AA
      CALL MXVDIR(N,-DIS,XO,S,XO)
      DIS=1.0D 0/(B*DIS)
      CALL MXDSMU(N,H,DIS,XO)
      CALL MXDSMU(N,H,-DIS,S)
      GO TO 21
   20 CONTINUE
*
*     GENERAL UPDATE
*
      DEN=PAR+POM*AA
      DIS=POM/DEN
      CALL MXDSMU(N,H,(PAR*DIS-1.0D 0)/A,S)
      CALL MXVDIR(N,-DIS,S,XO,S)
      CALL MXDSMU(N,H,DEN/B,S)
   21 CONTINUE
      IF (GAM.EQ.1.0D 0) GO TO 22
*
*     SCALING
*
      CALL MXDSMS(N,H,GAM)
   22 CONTINUE
      ITERH=0
      RETURN
      END
* SUBROUTINE PUDBM2                ALL SYSTEMS                92/12/01
C PORTABILITY : ALL SYSTEMS
C 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX.
*
* PARAMETERS :
*  RU  H(M)  POSITIVE DEFINITE APPROXIMATION OF THE HESSIAN
*         MATRIX.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  GO(NF)  GRADIENTS DIFFERENCE.
*
* COMMON DATA :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  M  NUMBER OF NONZERO ELEMENTS OF THE MATRIX.
*  II  MET  METHOD SELECTION. MET=1-BFGS UPDATE. MET=2-DFP UPDATE.
*         MET=3-HOSHINO UPDATE.
*  II  MET1  SELECTION OF SELF SCALING.  MET1=1-SELF SCALING SUPPRESSED.
*         MET1=2-INITIAL SELF SCALING.
*  II  MET2  SELECTION OF THE LINE SEARCH MODEL. MET2=1-QUADRATIC MODEL.
*         MET2=2 USE OF TAYLOR EXPANSION.
*  II  MET3  METHOD CORRECTION. MET3=1-NO CORRECTION.
*         MET3=2-POWELL'S CORRECTION.
*  II  ITERD  TERMINATION INDICATOR. ITERD<0-BAD DECOMPOSITION.
*         ITERD=0-DESCENT DIRECTION. ITERD=1-NEWTON LIKE STEP.
*         ITERD=2-INEXACT NEWTON LIKE STEP. ITERD=3-BOUNDARY STEP.
*         ITERD=4-DIRECTION WITH THE NEGATIVE CURVATURE.
*         ITERD=5-MARQUARDT STEP.
*  IO  ITERH  TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
*         ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
*  II  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*         IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=2-BUNCH-PARLETT
*         DECOMPOSITION. IDECF=3-INVERSION.
*  II  ITRAN  TRANSFORMATION INDICATOR. ITRAN=0 OR ITRAN=1 IF
*         TRANSFORMATION IS NOT OR IS USED.
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RI  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  TO  TUXX  TEXT INFORMATION ON THE CORRECTION USED.
*
* SUBPROGRAMS USED :
*  S   MXDSMM  MATRIX-VECTOR PRODUCT.
*  S   MXDSMU  CORRECTION OF A DENSE SYMMETRIC MATRIX.
*  S   MXDSMS  SCALING OF A DENSE SYMMETRIC MATRIX.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF VECTORS.
*  S   MXVNEG  COPYING OF A VECTOR WITH THE CHANGE OF THE SIGN.
*  S   MXVSCL  SCALING OF A VECTOR.
*  S   UOU1D1  PRINT OF ENTRY TO VARIABLE METRIC UPDATE.
*  S   UOU1D2  PRINT OF EXIT FROM VARIABLE METRIC UPDATE.
*
* METHOD :
* BASIC VARIABLE METRIC METHODS.
*
      SUBROUTINE PUDBM2(NF,N,H,HH,S,XO,GO,SO,FO,PAR,MET1,MET3,IDECF,
     & ITERH)
      INTEGER NF,N,MET1,MET3,IDECF,ITERH
      DOUBLE PRECISION H(NF*(NF+1)/2),HH(NF*(NF+1)/2),S(NF),XO(NF),
     & GO(NF),SO(NF),FO,PAR
      DOUBLE PRECISION DEN,A,B,C,GAM,POM,MXVDOT
      LOGICAL L1
      DOUBLE PRECISION CON
      PARAMETER (CON=1.0D-8)
      IF (IDECF.NE.0) THEN
      ITERH=-1
      RETURN
      ENDIF
      L1=MET1.GE.2
C
C     DETERMINATION OF THE PARAMETERS B, C
C
      CALL MXDSMM(N,H,XO,S)
      CALL MXVDIF(N,GO,S,SO)
      IF (MET3.EQ.2) CALL MXVSCL(N,1.0D0/SQRT(FO),SO,SO)
      B=MXVDOT(N,XO,SO)
      IF (B.LE.0.0D0) L1=.FALSE.
      A=0.0D0
      CALL MXDSMM(N,HH,XO,S)
      C=MXVDOT(N,XO,S)
      IF (C.LE.0.0D0) L1=.FALSE.
      IF (L1) THEN
C
C     DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
C
      GAM=C/B
      ELSE
      GAM=1.0D0
      ENDIF
      PAR=GAM
C
C     RANK ONE UPDATE
C
      DEN=PAR*B-C
      IF (ABS(DEN).LE.CON*MAX(CON,ABS(PAR*B),ABS(C))) THEN
      IF (B.GT.0.0D0.AND.C.GT.0.0D0) GO TO 1
      ITERH=4
      RETURN
      ENDIF
      POM=PAR*B/DEN
      CALL MXVDIR(N,-PAR,SO,S,S)
      CALL MXDSMU(N,HH,1.0D0/DEN,S)
      GO TO 5
C
C     BFGS UPDATE
C
    1 POM=0.0D0
      CALL MXDSMU(N,HH,PAR/B,SO)
      IF (C.GT.0.0D0) CALL MXDSMU(N,HH,-1.0D0/C,S)
      GO TO 5
    5 ITERH=0
      IF (GAM.NE.1.0D0) THEN
      CALL MXDSMS(N,HH,1.0D0/GAM)
      ENDIF
      RETURN
      END
* SUBROUTINE PUDBQ1                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* BROYDEN GOOD UPDATE OF A RECTANGULAR MATRIX AFTER THE QR
* DECOMPOSITION.
*
* PARAMETERS :
*  II  N  NUMBER OF VARIABLES.
*  II  NA  NUMBER OF EQUATIONS.
*  RU  H(N*(N+1)/2)  UPDATED UPPER TRIANGULAR MATRIX.
*  RI  ETA2  PARAMETER WHICH CONTROLS A NONSINGULARITY
*  RU  AG(N*NA)  UPDATED RECTANGULAR MATRIX.
*  RA  S(N)  AUXILIARY VECTOR.
*  RI  XO(N)  VECTOR OF VARIABLES DIFFERENCE.
*  RI  AFO(NA)  RIGHT HAND SIDES DIFFERENCE.
*  II  MET  VARIABLE METRIC UPDATE.
*  IO  ITERH  UPDATE INDICATOR. ITERH=0-SUCCESSFUL UPDATE.
*         ITERH>0-UNSUCCESSFUL UPDATE.
*  IU  IDECA  DECOMPOSITION INDICATOR.
*  II  NDECA  NUMBER OF DECOMPOSITIONS.
*
* METHOD :
* VARIOUS VARIABLE METRIC UPDATES INCLUDING BFGS UPDATE.
*
      SUBROUTINE PUDBQ1(N,NA,H,ETA2,AG,S,XO,AFO,MET,ITERH,IDECA,
     & NDECA)
      INTEGER N,NA,MET,INF,ITERH,IDECA,NDECA
      DOUBLE PRECISION H(*),ETA2,AG(*),S(*),XO(*),AFO(*)
      DOUBLE PRECISION DEN,MXVDOT
      IF (MET.LE.0) RETURN
      IF (IDECA.EQ.0) THEN
*
*     QR DECOMPOSITION
*
      DEN=ETA2
      CALL MXDRQF(N,NA,AG,H)
      CALL MXDPRC(N,H,INF,DEN)
      NDECA=NDECA+1
      IDECA=1
      ELSE IF (IDECA.NE.1) THEN
      ITERH=-1
      GO TO 22
      ENDIF
*
*     THE GOOD BROYDEN UPDATE
*
      DEN=MXVDOT(N,XO,XO)
      IF (DEN.LE.0.0D 0)  THEN
      ITERH=2
      GO TO 22
      ENDIF
      CALL MXVCOP(N,XO,S)
   21 CONTINUE
      CALL MXVNEG(N,XO,XO)
      CALL MXDPRM(N,H,XO,1)
      CALL MXDRMD(N,NA,AG,XO,1.0D 0,AFO,AFO)
      CALL MXDRQU(N,NA,AG,H,1.0D 0/DEN,AFO,S,XO,INF)
      ITERH=0
   22 CONTINUE
      RETURN
      END
* SUBROUTINE PUDFM1                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX.
*
* PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RU  B(M)  POSITIVE DEFINITE APPROXIMATION OF THE HESSIAN MATRIX.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  GO(NF)  GRADIENTS DIFFERENCE.
*  RI  F  CURRENT VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RI  ETA5  TOLERANCE FOR A HYBRID METHOD.
*  RI  ETA9  MAXIMUM FOR REAL NUMBERS.
*  II  IPOM1  METHOD INDICATOR.
*  IO  IPOM2  INDICATOR FOR SCALING.
*  II  MET  METHOD SELECTION. MET=0-NO UPDATE. MET=1-BFGS UPDATE.
*  II  MET1  SELECTION OF SELF SCALING.  MET1=1-SELF SCALING SUPPRESSED.
*         MET1=2 SELF SCALING IN THE FIRST ITERATION AFTER RESTART.
*         MET1=3-SELF SCALING IN EACH ITERATION.
*
* COMMON DATA :
*  IO  ITERH  TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
*         ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
*  II  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*  TO  TUXX  TEXT INFORMATION ON THE CORRECTION USED.
*
* SUBPROGRAMS USED :
*  S   MXDSMM  MATRIX-VECTOR PRODUCT.
*  S   MXDSMU  CORRECTION OF A DENSE SYMMETRIC MATRIX.
*  S   MXDSMS  SCALING OF A DENSE SYMMETRIC MATRIX.
*  RF  MXVDOT  DOT PRODUCT OF VECTORS.
*  S   UOERR1  ERROR MESAGES.
*  S   UOU1D1  PRINT OF ENTRY TO VARIABLE METRIC UPDATE.
*  S   UOU1D2  PRINT OF EXIT FROM VARIABLE METRIC UPDATE.
*
* METHOD :
*  FLETCHER'S COMBINATION OF THE GAUSS-NEWTON AND THE BFGS METHODS.
*
      SUBROUTINE PUDFM1(N,B,S,XO,GO,F,FO,ETA5,IPOM1,IPOM2,MET1,IDECF,
     & ITERH)
      INTEGER N,IPOM1,IPOM2,MET1,IDECF,ITERH
      DOUBLE PRECISION B(N*(N+1)/2),S(N),XO(N),GO(N),F,FO,ETA5
      DOUBLE PRECISION MXVDOT
      DOUBLE PRECISION AB,BB,CB,GAM,PAR
      LOGICAL L1
      IF (IDECF.NE.0) THEN
      ITERH=-1
      RETURN
      ENDIF
      PAR=0.0D0
C
C     DETERMINATION OF THE PARAMETERS A,B,C
C
      BB=MXVDOT(N,XO,GO)
      IF (BB.LE.0.0D0)  THEN
      ITERH=2
      IPOM1=0
      RETURN
      ENDIF
      AB=0.0D0
      CALL MXDSMM(N,B,XO,S)
      CB=MXVDOT(N,XO,S)
      IF (CB.LE.0.0D0) THEN
      ITERH=3
      RETURN
      ENDIF
      L1=MET1.EQ.4.OR.MET1.EQ.3.AND.IPOM2.GE.1.OR.MET1.EQ.2
     & .AND.IPOM2.EQ.1
      IF (FO-F.GE.ETA5*FO) THEN
      IPOM1=0
      ELSE
      IPOM1=1
      ENDIF
      IF (L1)  THEN
C
C     DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
C
      GAM=CB/BB
      ELSE
      GAM=1.0D0
      ENDIF
C
C     BFGS UPDATE
C
      CALL MXDSMU(N,B,GAM/BB,GO)
      CALL MXDSMU(N,B,-1.0D0/CB,S)
      ITERH=0
      IPOM2=0
      IF (L1) THEN
      CALL MXDSMS(N,B,1.0D0/GAM)
      ENDIF
      RETURN
      END
* SUBROUTINE PUDRV1                ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DRIVER FOR HYBRID QUASI-NEWTON UPDATES.
*
* PARAMETERS:
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  FO  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RI  F  CURRENT VALUE OF THE OBJECTIVE FUNCTION.
*  RI  PO  PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
*  II  IPOM1  UPDATE SELECTION.
*  II  IPOM2  METHOD SELECTION.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  IREST  RESTART SPECIFICATION. IF IREST=0 DOES NOT HOLD THEN A
*         RESTART IS PERFORMED.
*
      SUBROUTINE PUDRV1(R,FO,F,PO,IPOM1,IPOM2,NRED,IREST)
      INTEGER IPOM1,IPOM2,NRED,IREST
      DOUBLE PRECISION R,FO,F,PO
      DOUBLE PRECISION POM
      DOUBLE PRECISION CON1,CON2
      PARAMETER(CON1=1.0D-1,CON2=1.0D-2)
      POM=(FO-F)/FO
      GO TO (1,2,3,4) IPOM2
    1 IREST=1
      IF (NRED.LE.0) THEN
      IPOM1=2
      IREST=0
      ELSE
      IPOM1=0
      ENDIF
      GO TO 5
    2 IREST=1
      IF(POM.GE.CON2) THEN
      IPOM1=0
      ELSE IF (F-FO.LE.R*PO) THEN
      IPOM1=0
      ELSE
      IPOM1=1
      IREST=0
      ENDIF
      GO TO 5
    3 IREST=1
      IF (NRED.LE.0) THEN
      IF (IPOM1.NE.1) THEN
      IPOM1=2
      IREST=0
      ELSE
      IPOM1=0
      ENDIF
      ELSE IF (POM.GE.CON2) THEN
      IPOM1=0
      ELSE IF (IPOM1.NE.2) THEN
      IPOM1=1
      IREST=0
      ELSE
      IPOM1=0
      ENDIF
      GO TO 5
    4 IREST=1
      IPOM1=0
    5 CONTINUE
      RETURN
      END
* SUBROUTINE PUDSD2                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* INITIATION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
*
* PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RU  H(N*(N+1)/2)  POSITIVE DEFINITE APPROXIMATION OF THE HESSIAN
*         MATRIX
*  RU  B(N*(N+1)/2)  POSITIVE DEFINITE APPROXIMATION OF THE HESSIAN
*         MATRIX
*  RI  F  CURRENT VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RI  ETA5  TOLERANCE FOR A HYBRID METHOD.
*  II  MET3  TYPE OF STRUCTURED UPDATE. MET3=1-STANDARD STRUCTURED
*         UPDATE. MET3=2-TOTALLY STRUCTURED UPDATE.
*
* COMMON DATA :
*  RU  RAN  RANDOM NUMBER.
*  II  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*  II  IREST  RESTART SPECIFICATION. IF IREST=0 DOES NOT HOLD THEN A
*         RESTART IS PERFORMED.
*  II  IDIR INDICATOR OF DIRECTION DETERMINATION. IDIR=0-BASIC
*         DETERMINATION. IDIR=1-DETERMINATION AFTER STEPSIZE
*         REDUCTION. IDIR=2-DETERMINATION AFTER STEPSIZE EXPANSION.
*  IU  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  TO  TUXX   TEXT INFORMATION ON THE RESTART USED.
*
* SUBPROGRAMS USED :
*  S   MXDSMI  GENERATION OF THE UNIT MATRIX.
*  S   MXDSDO  INITIATION OF A DIAGONAL MATRIX.
*  S   MXDSMA  DENSE SYMMETRIC MATRIX AUGMENTED
*              BY THE SCALED DENSE SYMMETRIC MATRIX.
*  S   MXDSMO  A SCALAR IS SET TO ALL ELEMENTS
*              OF A DENSE SYMMETRIC MATRIX.
*  S   MXDSMS  SCALING OF A DENSE SYMMETRIC MATRIX.
*  S   UYSET3  DEFINITION OF THE RESTART VARIABLES.
*
      SUBROUTINE PUDSD2(N,H,B,F,FO,ETA5,MET3,LD,IDIR,IDECF,IREST,IND)
      INTEGER N,MET3,LD,IDIR,IDECF,IREST,IND
      DOUBLE PRECISION H(N*(N+1)/2),B(N*(N+1)/2),F,FO,ETA5
      INTEGER IUDSD
      SAVE IUDSD
      IND=0
      IF (IREST.LT.0) THEN
      CALL MXDSMI(N,B)
      IF (F.LT.1.0D0) CALL MXDSMS(N,B,SQRT(F))
      IUDSD=1
      ELSE IF (IREST.EQ.0) THEN
      IF (IDIR.LE.0) THEN
      IF (FO-F.LE.ETA5*FO) THEN
      IF (MET3.EQ.2) THEN
      CALL MXDSMA(N,SQRT(F),B,H,H)
      ELSE
      CALL MXDSMA(N,1.0D0,B,H,H)
      ENDIF
      LD=MIN(LD,1)
      IUDSD=0
      ENDIF
      ENDIF
      ELSEIF(IUDSD.EQ.0) THEN
      IF (MET3.EQ.2) THEN
      CALL MXDSMA(N,-SQRT(F),B,H,H)
      ELSE
      CALL MXDSMA(N,-1.0D0,B,H,H)
      ENDIF
      CALL MXDSMI(N,B)
      IF (F.LT.1.0D0) CALL MXDSMS(N,B,SQRT(F))
      IUDSD=1
      ELSE
      CALL MXDSMI(N,H)
      LD=MIN(LD,1)
      IUDSD=1
      IND=1
      ENDIF
      IDECF=0
      RETURN
      END
* SUBROUTINE PUDSD3                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* INITIATION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
*
* PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RU  H(N*(N+1)/2)  FACTORIZATION H=L*D*TRANS(L) OF A POSITIVE
*         SEMIDEFINITE HESSIAN MATRIX.
*  RU  B(N*(N+1)/2)  FACTORIZATION B=L*D*TRANS(L) OF A POSITIVE
*         DEFINITE APPROXIMATION OF THE HESSIAN MATRIX.
*  IU  IPOM1  METHOD INDICATOR.
*  IU  IPOM2  INDICATOR FOR SCALING.
*
* COMMON DATA :
*  RU  RAN  RANDOM NUMBER.
*  II  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*  II  IREST  RESTART SPECIFICATION. IF IREST=0 DOES NOT HOLD THEN A
*         RESTART IS PERFORMED.
*  II  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP.
*  IU  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*
* SUBPROGRAMS USED :
*  S   MXDSMI  GENERATION OF THE UNIT MATRIX.
*  S   MXDSDO  INITIATION OF A DIAGONAL MATRIX.
*  S   UUDSMC  COPYING OF A DENSE SYMMETRIC MATRIX.
*  RF  UNRAN1  RANDOM NUMBER GENERATOR.
*  S   UYSET3  DEFINITION OF THE RESTART VARIABLES.
*
      SUBROUTINE PUDSD3(N,H,B,IPOM1,IPOM2,LD,IDECF,ITERS,IREST,IND)
      INTEGER N,IPOM1,IPOM2,LD,IDECF,ITERS,IREST,IND
      DOUBLE PRECISION H(N*(N+1)/2),B(N*(N+1)/2)
      INTEGER KDECF
      SAVE KDECF
      IND=0
      IF (IREST.EQ.0.OR.IREST.GT.0.AND.IPOM1.EQ.0) THEN
      ELSE
      CALL MXDSMI(N,B)
      IPOM2= 1
      KDECF=-1
      ENDIF
      IF (IPOM1.EQ.1) THEN
      IF (ITERS.GT.0.OR.IREST.GT.0) THEN
      CALL MXDSMC(N,B,H)
      LD=MIN(LD,1)
      IF (IPOM2.EQ.1) IDECF=KDECF
      IF (IREST.GT.0) IND=1
      ENDIF
      ELSE
      IF (IREST.LE.0) THEN
      ELSE
      IPOM1= 1
      CALL MXDSMC(N,B,H)
      LD=MIN(LD,1)
      IF (IPOM2.EQ.1) IDECF=KDECF
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PYADB4             ALL SYSTEMS                   98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* NEW LINEAR CONSTRAINTS OR NEW SIMPLE BOUNDS ARE ADDED TO THE ACTIVE
* SET. GILL-MURRAY FACTORIZATION OF THE TRANSFORMED HESSIAN MATRIX
* APPROXIMATION IS UPDATED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  IU  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  CF(NC)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RI  CFD(NC) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT FUNCTIONS.
*  IU  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RU  H(NF*(NF+1)/2)  GILL-MURRAY FACTORIZATION OF THE TRANSFORMED
*         HESSIAN MATRIX APPROXIMATION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  II  KBF  TYPE OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. KBF=1-ONE
*         SIDED SIMPLE BOUNDS. KBF=2-TWO SIDED SIMPLE BOUNDS.
*  II  KBC  TYPE OF CONSTRAINTS. KBC=0-NO CONSTRAINTS. KBC=1-CONSTRAINTS
*         WITH ONE SIDED BOUNDS. KBC=2-CONSTRAINTS WITH TWO SIDED
*         BOUNDS.
*  IU  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IO  IER  ERROR INDICATOR.
*  IO  ITERM  TERMINATION INDICATOR.
*
* COMMON DATA :
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*
* SUBPROGRAMS USED :
*  S   PLADB4  ADDITION OF A NEW ACTIVE CONSTRAINT.
*  S   PLNEWS  IDENTIFICATION OF ACTIVE UPPER BOUNDS.
*  S   PLNEWL  IDENTIFICATION OF ACTIVE LINEAR CONSTRAINRS.
*  S   PLDIRL  NEW VALUES OF CONSTRAINT FUNCTIONS.
*  S   MXVIND  CHANGE OF THE INTEGER VECTOR FOR CONSTRAINT ADDITION.
*
      SUBROUTINE PYADB4(NF,N,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,
     & CR,CZ,H,S,R,EPS7,EPS9,GMAX,UMAX,KBF,KBC,INEW,IER,ITERM)
      INTEGER NF,N,NC,IX(*),IC(*),ICA(*),KBF,KBC,INEW,IER,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*),
     & CG(*),CR(*),CZ(*),H(*),S(*),R,EPS7,EPS9,GMAX,UMAX
      INTEGER I,J,K,L,IJ,IK,KC,KJ,KK,LL
      DOUBLE PRECISION DEN,TEMP
      INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      IF (KBC.GT.0) THEN
      IF (R.NE.0.0D 0) CALL PLDIRL(NC,CF,CFD,IC,R,KBC)
      IF (INEW.NE.0) THEN
      IF (KBF.GT.0) THEN
      DO 1 I=1,NF
      INEW=0
      CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW)
      CALL PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,9,INEW,
     & NADD,IER)
      CALL MXVIND(IX,I,IER)
      IF (IER.LT.0) THEN
      ITERM=-15
      RETURN
      ENDIF
    1 CONTINUE
      ENDIF
      DO 2 KC=1,NC
      INEW=0
      CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      CALL PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,9,INEW,
     & NADD,IER)
      CALL MXVIND(IC,KC,IER)
      IF (IER.LT.0) THEN
      ITERM=-15
      RETURN
      ENDIF
    2 CONTINUE
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      K=0
      DO 7 L=1,NF
      IF (IX(L).GE.0) K=K+1
      INEW=0
      CALL PLNEWS(X,IX,XL,XU,EPS9,L,INEW)
      IF (INEW.NE.0) THEN
      IX(L)=10-IX(L)
      KK=K*(K-1)/2
      DEN=H(KK+K)
      IF (DEN.NE.0.0D 0) THEN
      IJ=0
      KJ=KK
      DO 4 J=1,N
      IF (J.LE.K) THEN
      KJ=KJ+1
      ELSE
      KJ=KJ+J-1
      ENDIF
      IF (J.NE.K) TEMP=H(KJ)/DEN
      IK=KK
      DO 3 I=1,J
      IF (I.LE.K) THEN
      IK=IK+1
      ELSE
      IK=IK+I-1
      ENDIF
      IJ=IJ+1
      IF (I.NE.K.AND.J.NE.K) H(IJ)=H(IJ)+TEMP*H(IK)
    3 CONTINUE
    4 CONTINUE
      ENDIF
      LL=KK+K
      DO 6 I=K+1,N
      DO 5 J=1,I
      LL=LL+1
      IF (J.NE.K) THEN
      KK=KK+1
      H(KK)=H(LL)
      ENDIF
    5 CONTINUE
    6 CONTINUE
      N=N-1
      ENDIF
    7 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PYFUT1                ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TERMINATION CRITERIA AND TEST ON RESTART.
*
* PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RI  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RI  UMAX  MAXIMUN ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  RO  GMAX  NORM OF THE TRANSFORMED GRADIENT.
*  RI  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*  RI  TOLX  LOWER BOUND FOR STEPLENGTH.
*  RI  TOLF  LOWER BOUND FOR FUNCTION DECREASE.
*  RI  TOLB  LOWER BOUND FOR FUNCTION VALUE.
*  RI  TOLG  LOWER BOUND FOR GRADIENT.
*  II  KD  DEGREE OF REQUIRED DERIVATIVES.
*  IU  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER RESTART.
*  II  MIT  MAXIMUM NUMBER OF ITERATIONS.
*  IU  NFV  ACTUAL NUMBER OF COMPUTED FUNCTION VALUES.
*  II  MFV  MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES.
*  IU  NFG  ACTUAL NUMBER OF COMPUTED GRADIENT VALUES.
*  II  MFG  MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES.
*  IU  NTESX  ACTUAL NUMBER OF TESTS ON STEPLENGTH.
*  II  MTESX  MAXIMUM NUMBER OF TESTS ON STEPLENGTH.
*  IU  NTESF  ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE.
*  II  MTESF  MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE.
*  II  ITES  SYSTEM VARIBLE WHICH SPECIFIES TERMINATION. IF ITES=0
*         THEN TERMINATION IS SUPPRESSED.
*  II  IRES1  RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
*         IRES1*N+IRES2 ITERATIONS.
*  II  IRES2  RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
*         IRES1*N+IRES2 ITERATIONS.
*  IU  IREST  RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*  IO  ITERM  TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX
*         UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF
*         UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER
*         BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND
*         FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF
*         ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF
*         COMPUTED FUNCTION VALUES.
*
      SUBROUTINE PYFUT1(N,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD,
     & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,IRES1,
     & IRES2,IREST,ITERS,ITERM)
      INTEGER N,KD,NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,
     & ITES,IRES1,IRES2,IREST,ITERS,ITERM
      DOUBLE PRECISION  F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLG,TOLB
      DOUBLE PRECISION  TEMP
      IF (ITERM.LT.0) RETURN
      IF (ITES .LE.0) GO TO 2
      IF (ITERS.EQ.0) GO TO 1
      IF (NIT.LE.0) FO=F+MIN(SQRT(ABS(F)),ABS(F)/1.0D 1)
      IF (F.LE.TOLB) THEN
      ITERM = 3
      RETURN
      ENDIF
      IF (KD.GT.0) THEN
      IF (GMAX.LE.TOLG.AND.UMAX.LE.TOLG) THEN
      ITERM = 4
      RETURN
      ENDIF
      ENDIF
      IF (NIT.LE.0) THEN
      NTESX = 0
      NTESF = 0
      ENDIF
      IF (DMAX.LE.TOLX) THEN
      ITERM = 1
      NTESX = NTESX + 1
      IF (NTESX .GE. MTESX) RETURN
      ELSE
      NTESX = 0
      ENDIF
      TEMP=ABS(FO-F)/MAX(ABS(F),1.0D 0)
      IF (TEMP.LE.TOLF) THEN
      ITERM = 2
      NTESF = NTESF+1
      IF (NTESF.GE.MTESF) RETURN
      ELSE
      NTESF = 0
      ENDIF
    1 IF (NIT.GE.MIT) THEN
      ITERM = 11
      RETURN
      ENDIF
      IF (NFV.GE.MFV) THEN
      ITERM = 12
      RETURN
      ENDIF
      IF (NFG.GE.MFG) THEN
      ITERM = 13
      RETURN
      ENDIF
    2 ITERM = 0
      IF (N.GT.0.AND.NIT-KIT.GE.IRES1*N+IRES2) THEN
      IREST=MAX(IREST,1)
      ENDIF
      NIT = NIT + 1
      RETURN
      END
* SUBROUTINE PYRMB1               ALL SYSTEMS                98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* OLD LINEAR CONSTRAINT OR AN OLD SIMPLE BOUND IS REMOVED FROM THE
* ACTIVE SET. TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION AND
* TRANSFORMED HESSIAN MATRIX APPROXIMATION ARE UPDATED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  IU  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  IU  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  GN(NF)  TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  H(NF*(NF+1)/2)  TRANSFORMED HESSIAN MATRIX APPROXIMATION.
*  RI  EPS8  TOLERANCE FOR CONSTRAINT TO BE REMOVED.
*  RI  UMAX  MAXIMUN ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  RI  GMAX  NORM OF THE TRANSFORMED GRADIENT.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  IOLD  INDEX OF THE REMOVED CONSTRAINT.
*  IA  KOLD  AUXILIARY VARIABLE.
*  IA  KREM  AUXILIARY VARIABLE.
*  IO  IER  ERROR INDICATOR.
*  IO  ITERM  TERMINATION INDICATOR.
*
* COMMON DATA :
*  IU  NREM  NUMBER OF CONSTRAINT DELETIONS.
*
* SUBPROGRAMS USED :
*  S   PLRMB0  CONSTRAINT DELETION.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PYRMB1(NF,N,IX,IC,ICA,CG,CR,CZ,G,GN,H,EPS8,UMAX,
     & GMAX,KBF,KBC,IOLD,KOLD,KREM,IER,ITERM)
      INTEGER NF,N,IX(*),IC(*),ICA(*),KBF,KBC,IOLD,KOLD,KREM,IER,
     $ ITERM
      DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),H(*),EPS8,UMAX,
     & GMAX
      INTEGER I,J,K,KC,L
      INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      IF (KBC.GT.0) THEN
      IF (UMAX.GT.EPS8*GMAX) THEN
      CALL PLRMB0(NF,N,ICA,CG,CR,CZ,G,GN,IOLD,KREM,NREM,IER)
      IF (IER.LT.0) THEN
      ITERM=-16
      ELSE IF (IER.GT.0) THEN
      IOLD=0
      ELSE
      K=N*(N-1)/2
      CALL MXVSET(N,0.0D 0,H(K+1))
      H(K+N)=1.0D 0
      KC=ICA(NF-N+1)
      IF (KC.GT.0) THEN
      IC(KC)=-IC(KC)
      ELSE
      K=-KC
      IX(K)=-IX(K)
      ENDIF
      ENDIF
      ELSE
      IOLD=0
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      IF (UMAX.GT.EPS8*GMAX)  THEN
      IX(IOLD)=MIN(ABS(IX(IOLD)),3)
      DO 1 I=N,KOLD,-1
      GN(I+1)=GN(I)
    1 CONTINUE
      GN(KOLD)=G(IOLD)
      N=N+1
      K=N*(N-1)/2
      L=K+N
      DO 3 I=N,KOLD,-1
      DO 2 J=I,1,-1
      IF (I.NE.KOLD.AND.J.NE.KOLD) THEN
      H(L)=H(K)
      K=K-1
      L=L-1
      ELSE IF (I.EQ.KOLD.AND.J.EQ.KOLD) THEN
      H(L)=1.0D 0
      L=L-1
      ELSE
      H(L)=0.0D 0
      L=L-1
      ENDIF
    2 CONTINUE
    3 CONTINUE
      ELSE
      IOLD=0
      KOLD=0
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRBD             ALL SYSTEMS                   98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
* AND TRANSFORMED. TEST VALUE DMAX IS DETERMINED.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  GO(NF)  GRADIENTS DIFFERENCE.
*  RI  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM CURRENT
*         REDUCED SUBSPACE.
*  RU  SN(NF)  TRANSFORMED DIRECTION VECTOR.
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RU  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RU  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RU  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*
* SUBPROGRAMS USED :
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY TRANSPOSE OF A DENSE
*         RECTANGULAR MATRIX.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*  S   MXVSCL  SCALING OF A VECTOR.
*
      SUBROUTINE PYTRBD(NF,N,X,IX,XO,G,GO,CZ,SN,R,F,FO,P,PO,DMAX,ITERS,
     & KBF,KBC)
      INTEGER NF,N,IX(*),ITERS,KBF,KBC
      DOUBLE PRECISION X(*),XO(*),G(*),GO(*),CZ(*),SN(*),R,F,FO,P,PO,
     & DMAX
      INTEGER I,K
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NF,G,GO,GO)
      PO=R*PO
      P=R*P
      ELSE
      F = FO
      P = PO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NF,G,GO)
      ENDIF
      DMAX = 0.0D 0
      IF (KBC.GT.0) THEN
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
    1 CONTINUE
      IF (N.GT.0) THEN
      CALL MXVSCL(N,R,SN,XO)
      CALL MXVCOP(NF,GO,SN)
      CALL MXDRMM(NF,N,CZ,SN,GO)
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      K=0
      DO 2 I=1,NF
      IF (IX(I).LT.0) GO TO 2
      K=K+1
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
      XO(K)=XO(I)
      GO(K)=GO(I)
    2 CONTINUE
      ELSE
      DO 3 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
    3 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRBG               ALL SYSTEMS                98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED.
* TEST VALUES GMAX AND UMAX ARE COMPUTED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GN(NF)  TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  RO  GMAX  NORM OF THE TRANSFORMED GRADIENT.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  IOLD  INDEX OF THE REMOVED CONSTRAINT.
*  IA  KOLD  AUXILIARY VARIABLE.
*
* SUBPROGRAMS USED :
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDPRB  BACK SUBSTITUTION.
*  S   MXVCOP  COPYING OF A VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*  RF  MXVMAX  L-INFINITY NORM OF A VECTOR.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*
      SUBROUTINE PYTRBG(NF,N,IX,IC,ICA,CG,CR,CZ,G,GN,UMAX,GMAX,
     & KBF,KBC,IOLD,KOLD)
      INTEGER NF,N,IX(*),IC(*),ICA(*),KBF,KBC,IOLD,KOLD
      DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),UMAX,GMAX
      DOUBLE PRECISION TEMP,MXVMAX,MXVDOT
      INTEGER NCA,NCZ,I,J,K,KC
      IOLD=0
      KOLD=0
      UMAX=0.0D 0
      GMAX=0.0D 0
      IF (KBC.GT.0) THEN
      IF (NF.GT.N) THEN
      NCA=NF-N
      NCZ=N*NF
      CALL MXVCOP(NF,G,GN)
      DO 1 J=1,NCA
      K=ICA(J)
      IF (K.GT.0) THEN
      CZ(NCZ+J)=MXVDOT(NF,CG((K-1)*NF+1),GN)
      ELSE
      I=-K
      CZ(NCZ+J)=GN(I)
      ENDIF
    1 CONTINUE
      CALL MXDPRB(NCA,CR,CZ(NCZ+1),0)
      DO 2 J=1,NCA
      TEMP=CZ(NCZ+J)
      KC=ICA(J)
      IF (KC.GT.0) THEN
      K=IC(KC)
      ELSE
      I=-KC
      K=IX(I)
      ENDIF
      IF (K.LE.-5) THEN
      ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN
      ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN
      ELSE
      IOLD=J
      UMAX=ABS(TEMP)
      ENDIF
    2 CONTINUE
      ENDIF
      IF (N.GT.0) THEN
      CALL MXDRMM(NF,N,CZ,G,GN)
      GMAX=MXVMAX(N,GN)
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      J=0
      IOLD=0
      KOLD=0
      DO 3 I=1,NF
      TEMP=G(I)
      K=IX(I)
      IF (K.GE.0) THEN
      J=J+1
      GN(J)=TEMP
      GMAX=MAX(GMAX,ABS(TEMP))
      ELSE IF (K.LE.-5) THEN
      ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN
      ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN
      ELSE
      IOLD=I
      KOLD=J+1
      UMAX=ABS(TEMP)
      ENDIF
    3 CONTINUE
      N=J
      ELSE
      DO 4 I=1,NF
      TEMP=G(I)
      GMAX=MAX(GMAX,ABS(TEMP))
    4 CONTINUE
      N=NF
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRBH               ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* HESSIAN MATRIX OF THE OBJECTIVE FUNCTION OR ITS APPROXIMATION IS
* SCALED AND REDUCED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RI  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  H(NF*(NF+1)/2)  HESSIAN MATRIX OR ITS APPROXIMATION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*
* SUBPROGRAMS USED :
*  S   MXDSMM  MATRIX VECTOR PRODUCT.
*  S   MXVCOP  COPYING OF A VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PYTRBH(NF,N,IX,CR,CZ,H,S,KBF,KBC,LD,ITERS)
      INTEGER NF,N,IX(*),KBF,KBC,LD,ITERS
      DOUBLE PRECISION CR(*),CZ(*),H(*),S(*)
      DOUBLE PRECISION MXVDOT
      INTEGER NCA,NCR,ICZ,JCZ,I,J,K,L
      IF (LD.NE.2.OR.ITERS.EQ.0) RETURN
      IF (KBC.GT.0) THEN
      IF (N.LE.0) RETURN
      NCA=NF-N
      NCR=NCA*(NCA+1)/2
      K=NCR
      JCZ=1
      DO 4 J=1,N
      CALL MXDSMM(NF,H,CZ(JCZ),S)
      ICZ=1
      DO 3 I=1,J
      K=K+1
      CR(K)=MXVDOT(NF,CZ(ICZ),S)
      ICZ=ICZ+NF
    3 CONTINUE
      JCZ=JCZ+NF
    4 CONTINUE
      CALL MXVCOP(N*(N+1)/2,CR(NCR+1),H)
      ELSE IF (KBF.GT.0) THEN
      K=0
      L=0
      DO 2 I=1,NF
      DO 1 J=1,I
      K=K+1
      IF(IX(I).LT.0.OR.IX(J).LT.0) GO TO 1
      L=L+1
      H(L)=H(K)
    1 CONTINUE
    2 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRBS               ALL SYSTEMS                98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED.
* VECTORS X,G AND VALUES F,P ARE SAVED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GO(NF)  SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS.
*  RO  CFD(NF)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  SN(NF)  TRANSFORMED DIRECTION VECTOR.
*  RO  S(NF)  DIRECTION VECTOR.
*  RO  RO  SAVED VALUE OF THE STEPSIZE PARAMETER.
*  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RU  FO  SAVED VALUE OF THE OBJECTIVE FUNCTION.
*  RI  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  PO  SAVED VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RU  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  IO  KREM  INDICATION OF LINEARLY DEPENDENT GRADIENTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE FUNCTION.
*
* SUBPROGRAMS USED :
*  S   PLMAXS  DETERMINATION OF THE MAXIMUM STEPSIZE USING SIMPLE
*         BOUNDS.
*  S   PLMAXL  DETERMINATION OF THE MAXIMUM STEPSIZE USING LINEAR
*         CONSTRAINTS.
*  S   MXDCMM  MATRIX VECTOR PRODUCT.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PYTRBS(NF,N,NC,X,IX,XO,XL,XU,G,GO,CF,CFD,IC,CL,CU,CG,
     & CZ,SN,S,RO,FP,FO,F,PO,P,RMAX,KBF,KBC,KREM,INEW)
      INTEGER NF,N,NC,IX(*),IC(*),KBF,KBC,KREM,INEW
      DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),G(*),GO(*),CF(*),CFD(*),
     & CL(*),CU(*),CG(*),CZ(*),SN(*),S(*),RO,FP,FO,F,PO,P,RMAX
      INTEGER I,K
      FP = FO
      RO = 0.0D 0
      FO = F
      PO = P
      CALL MXVCOP(NF,X,XO)
      CALL MXVCOP(NF,G,GO)
      IF (KBC.GT.0) THEN
      IF (N.GT.0) THEN
      CALL MXDCMM(NF,N,CZ,SN,S)
      INEW=0
      CALL PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,RMAX,KBC,KREM,INEW)
      CALL PLMAXS(NF,X,IX,XL,XU,S,RMAX,KBF,KREM,INEW)
      ELSE
      CALL MXVSET(NF,0.0D 0,S)
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      K=N+1
      DO 1 I=NF,1,-1
      IF (IX(I).LT.0) THEN
      S(I)=0.0D 0
      ELSE
      K=K-1
      S(I)=SN(K)
      ENDIF
    1 CONTINUE
      INEW=0
      CALL PLMAXS(NF,X,IX,XL,XU,S,RMAX,KBF,KREM,INEW)
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRFD             ALL SYSTEMS                   90/12/01
* PORTABILITY : ALL SYSTEMS
* 90/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* PREPARATION OF VARIABLE METRIC UPDATE.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RU  XO(NF)  SAVED VECTOR OF VARIABLES.
*  II  IAA(NF+1)  VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS.
*  RI  AG(NF*NA)  MATRIX WHOSE COLUMNS ARE GRADIENTS OF THE LINEAR
*          APPROXIMATED FUNCTIONS.
*  RI  AZ(NF+1)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  G(NF)  GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RU  GO(NF)  SAVED GRADIENT OF THE LAGRANGIAN FUNCTION.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IU  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  RU  R  VALUE OF THE STEPSIZE PARAMETER.
*  RU  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  SAVED VALUE OF THE OBJECTIVE FUNCTION.
*  RU  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RU  PO  SAVED VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  DMAX  RELATIVE STEPSIZE.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
*         LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
*         STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
*         ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
*         ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
*         ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
*         DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRFD(NF,NC,X,XO,IAA,AG,AZ,CG,G,GO,N,KD,LD,R,F,FO,P,
     +                  PO,DMAX,ITERS)
      DOUBLE PRECISION DMAX,F,FO,P,PO,R
      INTEGER          ITERS,KD,LD,N,NC,NF
      DOUBLE PRECISION AG(*),AZ(*),CG(*),G(*),GO(*),X(*),XO(*)
      INTEGER          IAA(*)
      INTEGER          I,J,L
      CALL MXVSET(NF,0.0D0,G)
      DO 10 J = 1,NF - N
          L = IAA(J)
          IF (L.GT.NC) THEN
              L = L - NC
              CALL MXVDIR(NF,-AZ(J),AG((L-1)*NF+1),G,G)
          ELSE IF (L.GT.0) THEN
              CALL MXVDIR(NF,-AZ(J),CG((L-1)*NF+1),G,G)
          ELSE
              L = -L
              G(L) = G(L) - AZ(J)
          END IF
   10 CONTINUE
      IF (ITERS.GT.0) THEN
          CALL MXVDIF(NF,X,XO,XO)
          CALL MXVDIF(NF,G,GO,GO)
          PO = R*PO
          P = R*P
      ELSE
          R = 0.0D0
          F = FO
          P = PO
          CALL MXVSAV(NF,X,XO)
          CALL MXVSAV(NF,G,GO)
          LD = KD
      END IF
      DMAX = 0.0D0
      DO 20 I = 1,NF
          DMAX = MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0))
   20 CONTINUE
      N = NF
      RETURN
      END
* SUBROUTINE PYTRND             ALL SYSTEMS                   91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX
* APPROXIMATION.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC NUMBER OF CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  XN(NF)  VECTOR OF SCALING FACTORS.
*  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RO  CZ(NF)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RO  CZS(NF)  SAVED VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  G(NF)  GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RI  GO(NF)  SAVED GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  CMAX  VALUE OF THE CONSTRAINT VIOLATION.
*  RO  CMAXO  SAVED VALUE OF THE CONSTRAINT VIOLATION.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*
* COMMON DATA :
*  II  NORMF  SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED.
*         NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY.
*         NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRND(NF,N,X,XO,ICA,CG,CZ,G,GO,R,F,FO,P,PO,CMAX,CMAXO,
     & DMAX,KD,LD,ITERS)
      INTEGER NF,N,KD,LD,ITERS
      INTEGER ICA(*)
      DOUBLE PRECISION X(*),XO(*),CG(*),CZ(*),G(*),GO(*),R,F,FO,P,PO,
     & CMAX,CMAXO,DMAX
      INTEGER I,J,L
      DO 2 J=1,NF-N
      L=ICA(J)
      IF (L.GT.0) THEN
      CALL MXVDIR(NF,-CZ(J),CG((L-1)*NF+1),G,G)
      ELSE
      L=-L
      G(L)=G(L)-CZ(J)
      ENDIF
    2 CONTINUE
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NF,G,GO,GO)
      PO=R*PO
      P=R*P
      ELSE
      F = FO
      P = PO
      CMAX = CMAXO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NF,G,GO)
      LD = KD
      ENDIF
      DMAX = 0.0D0
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0))
    1 CONTINUE
      N=NF
      RETURN
      END
* SUBROUTINE PYTRUD             ALL SYSTEMS                   98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
* AND SCALED. TEST VALUE DMAX IS DETERMINED.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  GO(NF)  GRADIENTS DIFFERENCE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*
* SUBPROGRAMS USED :
*  S   PYSET1  DEGREE DEFINITION OF THE COMPUTED DERIVATIVES.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRUD(NF,X,XO,G,GO,R,F,FO,P,PO,DMAX,KD,LD,ITERS)
      INTEGER NF,KD,LD,ITERS
      DOUBLE PRECISION X(*),XO(*),G(*),GO(*),R,F,FO,P,PO,DMAX
      INTEGER I
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NF,G,GO,GO)
      PO=R*PO
      P=R*P
      ELSE
      F = FO
      P = PO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NF,G,GO)
      LD=KD
      ENDIF
      DMAX = 0.0D 0
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PYTRUF             ALL SYSTEMS                   98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTORS OF VARIABLES DIFFERENCE AND RIGHT HAND SIDES DIFFERENCE ARE
* COMPUTED AND SCALED. TEST VALUE DMAX IS DETERMINED.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  NA NUMBER OF APPROXIMATED FUNCTIONS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  AF(NA)  VECTOR OF RIGHT HAND SIDES.
*  RI  AFO(NA)  VECTOR OF RIGHT HAND SIDES DIFFERENCE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*
* SUBPROGRAMS USED :
*  S   PYSET1  DEGREE DEFINITION OF THE COMPUTED DERIVATIVES.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRUF(NF,NA,X,XO,AF,AFO,R,F,FO,P,PO,DMAX,KD,LD,
     & ITERS)
      INTEGER NF,NA,KD,LD,ITERS
      DOUBLE PRECISION X(*),XO(*),AF(*),AFO(*),R,F,FO,P,PO,DMAX
      INTEGER I
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NA,AF,AFO,AFO)
      PO=R*PO
      P=R*P
      ELSE
      R = 0.0D 0
      F = FO
      P = PO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NA,AF,AFO)
      LD=KD
      ENDIF
      DMAX = 0.0D 0
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
    1 CONTINUE
      RETURN
      END
